An Historical Review of the Simplified Physical Fire Spread Model PhyFire: Model and Numerical Methods

Abstract

A historical review is conducted of PhyFire, a simplified physical forest fire spread model developed by the research group on Numerical Simulation and Scientific Computation (SINUMCC) at the University of Salamanca. The review ranges from the first version of the model to the current one now integrated into GIS, considering all the mathematical problems and numerical methods involved throughout its development: finite differences, mixed, classical and adaptive finite elements, data assimilation, sensitivity analysis, parameter adjustment, and parallel computation, among others. The simulation of processes as complex as forest fires involves a multidisciplinary effort that is constantly being enhanced, while posing interesting challenges from a mathematical, numerical, and computational perspective, without losing sight of the overriding aim of developing an efficient, effective, and useful simulation tool.

1. Introduction

The EU figure of 770,710 hectares burned in 2022 is the highest since records began in 2006. Spain leads the ranking with 299,436 hectares devastated by fire, more than triple the ten-year average, according to data provided by the European Forest Fire Information System (EFFIS) consulted at the end of September (see Figure 1 and [1]). In just one province, Zamora, more than 60,000 hectares have been burned in two of the largest fires thus far in Spain.

Not only were there major fires in Southern Europe in 2022, but also further north in countries that are not usually affected. It has been a year of all-time records, including heat waves and severe drought. The catastrophic figure of thousands of hectares devastated by fire in Europe alone shows that the problem is increasingly pressing under the growing threat of climate change.

This increase in the number and magnitude of wildfires also influences the accompanying smoke emissions. According to data from the CAMS Global Fire Assimilation System (GFAS), total wildfire emissions in the EU during the summer of 2022 are estimated to have amounted to 6.4 megatons of carbon, the highest level since 2007 [2].

This situation has not only affected Europe, as in the first half of 2022, extreme fires have also ravaged vast swathes of land across the world. In late July and early August, numerous forest fires broke out in the eastern and western regions of Eurasia. The large number of wildfires in 2021 in Western U.S. states made it the most-devastating year in this area so far. Their rate has nonetheless increased in the U.S. Pacific Northwest and Western Canada.

The above data provided only a very partial view of the problem of forest fires. Africa, where fires are widely used to clear grasslands and savanna, accounts for around seventy percent of the area burned globally. Many of the fires in South America involve the Amazon rainforest. The latest fire season in the Amazon region (August–September) has registered the highest number of fires in the last decade, although it is not clear whether this upward trend will continue. The bushfire season in Australia recorded the most-serious fires. For an overall view of the global wildfire situation, we recommend [3] and its references.

Although some wildfires occur naturally and help keep forest ecosystems healthy, an astonishing percentage are caused by human-related activities. Wildfires not only destroy lives, wildlife habitats, forests, and property, they also accelerate climate change.

Understanding and predicting the behavior of a system as complex as a wildfire is an undeniably useful tool for reducing its negative effects. Mathematical modeling and numerical simulation play a fundamental role, and they can ultimately assist in decision-making in prevention, fire-fighting, and further analysis. Improved remote sensing information technologies, advanced computational capabilities, and the development of communication technology have exponentially increased the efficiency and applicability of wildfire modeling. The quasi-empirical models combined with the software developed for extending the utility of these models to the landscape level has resulted in comprehensive and accurate tools for predicting wildfire spread, such as FlamMap-FARSITE [4,5] and Prometheus [6].

The rapid increase in computing power and advances in technology allow more complex models to become a real option, so research is focusing on physical-based models coupled with the atmosphere [7], for example HIGRAD/FIRETEC [8] or WFDS [9]. These models are computationally expensive, but there are also other more affordable atmosphere–wildfire models combining empirical or simplified quasi-physical fire models with atmospheric ones, such as WRF-SFIRE [10] or MesoNH-ForeFire [11].

The advances made in remote sensing information technologies (i.e., satellites, aircraft, and drones) have provided more and better data for feeding wildfire simulation models. A new type of model, referred to as data-driven or data-assimilation [12], is a result of these technologies. Their objective is to reduce the uncertainties in both model fidelity and input data by using real-time observations of wildfire dynamics. Artificial intelligence can also exploit the current ability to obtain increasingly more reliable data on wildfires for predicting their behavior [13]. Nevertheless, the use of wildfire spread modeling and the tools developed accordingly have been relatively limited in practice due to the high level of uncertainty in both the models and the available data for their initialization and parametrization. Reducing the level of uncertainty continues to be one of the unresolved challenges of wildfire modeling.

The model described in this paper, PhyFire, is framed within the physical-based models. It is a single-phase simplified two-dimensional model based on the fundamental physics of combustion and fire spread. It considers convection and radiation as the main heat transfer mechanisms, taking into account the heat lost by natural convection, the effect of the flame tilt caused by wind or slope, the influence of fuel moisture content and fuel type, and certain random effects such as fire-spotting. The resulting Partial Differential Equations (PDEs) are solved using efficient numerical and computational tools to obtain a software with levels of efficiency comparable to empirical models, including its own finite element toolbox and the use of parallel computing techniques.

Furthermore, PhyFire can be coupled with HDWind [14,15], which is a high-resolution wind field that adjusts a set of specific wind measurements to local factors within the simulation domain, such as the slope, roughness of the terrain, and surface temperature gradients. This wind model is based on an asymptotic approximation of the Navier–Stokes equations, deriving a mass-consistent vertical diffusion model capable of providing a 3D wind field by solving only 2D linear equations. This model was developed at the same time as the PhyFire model, by the same authors, with the initial aim of improving wind data for the fire spread model, but it was soon realized that the wind model posed very interesting mathematical and computational challenges and provided new applications [16].

PhyFire and HDWind have now been integrated into a Geographical Information System (GIS) [17], which is available for use in Spain through the URL: http://sinumcc.usal.es (accessed on 30 January 2023) [18].

This paper covers the PhyFire design process from its inception to the current GIS-integrated version, emphasizing its mathematical and numerical aspects and explaining why the modeling of such complex questions is a major source of interesting mathematical and numerical problems. Some of the results included here are unpublished, such as the flame temperature sub-model, the new predictor–corrector numerical scheme, and the use of new pre- and post-processing functionalities.

The paper is organized as follows. Section 2 describes the initial fire spread models in chronological order, through their main developments, focusing on mathematical and numerical aspects. Section 3 summarizes the current PhyFire model and describes in detail the numerical method used, as well as some original ways of improving the model efficiency. Section 4 discusses certain details of the work performed to integrate the model into a GIS. Section 5 covers a real wildfire simulated with the PhyFire-HDWind system. We end with the conclusions and suggestions for future work.

2. Previous Models

2.1. Conservation Laws

Combustion in a forest fire is a physical–chemical process that, as such, can be represented by conservation laws. These laws can be written in general as PDEs having the following form:

$$\frac{\partial }{\partial t}\left(\rho \varphi \right)+\nabla \left(\rho \varphi \mathbf{v}\right)-\nabla \left(K\nabla \varphi \right)=S_{\varphi }$$

(1)

where $\varphi$ is a generic magnitude, $\rho$ represents the density, $\mathbf{v}=(v_x,v_y,v_z)$ is the velocity, K is the diffusion coefficient, and $S_{\varphi }$ is the source term. Initially, we considered a three-dimensional problem, that is $\nabla =(\frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial z})$.

A first paramount simplification to achieve an efficient simulation model from the computational point of view is to reduce the three-dimensional problem to a two-dimensional one. Equation (1) describes the evolution of the corresponding magnitude in a three-dimensional domain, which represents a bed of fuel of thickness $\delta z$. We considered the following boundary conditions: homogeneous Dirichlet on the lateral boundary of the three-dimensional domain, homogeneous Neumann on the lower boundary $(z=0)$, and non-homogeneous Neumann on the upper boundary $(z=\delta z)$:

$$K\frac{\partial \varphi }{\partial \eta }{|}_{{\Gamma }_{z=\delta z}}=H_{\varphi }.$$

Denoting

$$\overline{\varphi }(x,y)={\int }_0^{\delta z}\varphi (x,y,z)\mathrm{d}z$$

and integrating Equation (1) with respect to the variable z, we can reduce the problem to two dimensions in space.

$$\frac{\partial }{\partial t}\left(\rho \overline{\varphi }\right)+\nabla \left(\rho \overline{\varphi }\mathbf{v}\right)-\nabla (K_{\varphi }\nabla \overline{\varphi })={\overline{S}}_{\varphi }+H_{\varphi },$$

(2)

where, now, $\nabla =(\frac{\partial }{\partial x},\frac{\partial }{\partial y})$ and $\mathbf{v}=(v_x,v_y)$.

We assumed that, $v_z=0$ away from the fire, the vertical variation of the corresponding magnitude is null near the fire, $K_{\varphi }$ is not dependent on z, and the thickness of the three-dimensional domain $\delta z$ is small enough to suppose that ${\overline{S}}_{\varphi }\approx S_{\overline{\varphi }}$.

Another aspect that adds complexity to an eventual fire spread simulation model is the fact that combustion in a forest fire is a positive feedback system, so we can distinguish two phases, an endothermic phase or solid phase, in which heat is used to release the volatile substances, and an exothermic phase or gas phase, in which the volatile substances mixed with oxygen react, producing more heat (see [19]). Bearing this in mind, a two-phase model can be proposed, with an energy conservation equation for each phase and mass conservation equations for solid fuel, gaseous fuel, and oxygen [20]. A two-phase model poses the drawback of correctly evaluating nonlinear terms. The correct evaluation of these terms would lead to introducing more elaborate models of turbulence, leading to much more complex systems, moving away from the goal of designing an efficient computational model.

Each of the models that we reviewed depends on the simplifying hypotheses raised in each case. All of them attempt to represent the main heat transfer mechanisms in a wildfire, namely radiation and convection [21]. The predominant mode of heat transfer will depend on the wind conditions, terrain slope, fuel type, and location relative to the fire.

2.2. First Model: Turbulent Diffusivity

The first simplified physical model for wildfire spread simulation proposed was a two-dimensional one-phase model considering turbulent flow, vertical heat lost, and a convective term representing the effect of wind [22]. The original equations are

$$\rho C\left(\frac{\partial }{\partial t}T+\mathbf{v}·\nabla T\right)-\nabla \left(K\nabla T\right)+H\left(T\right)(T-T_{\infty })=QMA{e}^{-\frac{E_A}{RT}},$$

(3)

$$\frac{\partial }{\partial t}M=-MA{e}^{-\frac{E_A}{RT}},$$

(4)

where T is the average value of temperature and M is the solid fuel load. K is the turbulent diffusivity. $H\left(T\right)=\overline{H}{(T-T_{\infty })}^{1/3}$ is the vertical convection heat transfer. $T_{\infty }$ is the ambient temperature. The reaction rate is given by the Arrhenius law. R is the universal gas constant. $E_A$ represents the activation energy. A is the pre-exponential or frequency factor of the reaction. We assumed that the activation energy for the gaseous phase is much lower than the one of the solid phase in order to propose a simplified one-phase chemical reaction. Q is the global heat production in the reaction.

The PDEs to be solved are a non-dimensional version of Equations (3) and (4), based on the following change of variables: the dimensionless temperature is $u=(T-T_{\infty })/(\overline{T}-T_{\infty })$, where $\overline{T}$ is a temperature up to which generalized combustion takes place; the dimensionless solid fuel is $c=M/M_0$, where $M_0$ the initial solid fuel load. Dimensionless time and space variables were also introduced, although, for the sake of simplicity, we kept the same notation. The non-dimensional equations are

$$\frac{\partial }{\partial t}u+\mathbf{w}·\nabla u-\nabla (\kappa \nabla u)+h\left(u\right)u=qcs\left(u\right),$$

(5)

$$\frac{\partial }{\partial t}c=-cs\left(u\right),$$

(6)

where $\mathbf{w}$ is the dimensionless wind velocity field, $\kappa$ is related to the turbulent diffusivity, $h\left(u\right)$ is the normalized vertical convection coefficient, q depends on the heat released, and $s\left(u\right)$ is written in terms of the previous exponential expression.

The non-dimensional equations were completed with the initial and boundary conditions. We propose the Dirichlet boundary conditions assuming the domain is big enough to ensure that the fire does not arrive at the boundary. The initial conditions represent the initial fuel load and fire source.

The first numerical method proposed was an implicit and upwind finite difference scheme for the energy Equation (5), assuming that the wind field velocity and diffusivity are constant along the domain. The matrix of the linear system obtained in each time step is block tridiagonal and strictly diagonally dominant, so the blockwise Gauss–Seidel method is convergent. The global Lipschitz condition that the source term satisfies ensures the stability of the scheme. An implicit Euler method was the scheme proposed for the fuel Equation (6).

The second numerical method proposed for the partial differential problem given by Equations (5) and (6) was an Adaptive Finite-Element Method (AFEM). The idea of using an adaptive method aims to reduce the operational cost by adapting the mesh to the numerical solution where more precision is necessary, the fire front, as this takes up a small part of the whole domain. This method combines refinement and derefinement techniques to generate in each time step a sequence of nested meshes, allowing an easy application of multigrid acceleration techniques. The design of this first model and its numerical resolution already pursued the aim of simulating the fire spread in times much shorter than real-time, through a simplified model and efficient numerical methods.

2.3. Second Model: Rosseland Approximation for Local Radiation

The second model is again a one-phase 2D model, but it addresses radiation as the dominant thermal transfer mechanism, without forgetting convection, which represents the effect of wind and slope. In [23], radiation was incorporated in the fire spread model by using a local radiation term using the Stefan–Boltzmann law and approximating the fourth power of the temperature by a Taylor expansion [24]. This is known as the Rosseland radiation model; it is valid when the medium is optically thick. It is a diffusion approximation for the radiation, and it introduces a nonlinearity in the diffusive term of the energy equation. Another modification compared to the first model is the simplification of the vertical heat loss term or natural convection. In this second model, the one-phase simplification leads to introducing a logic expression depending on the phase change temperature, which we can assimilate with the pyrolysis temperature. The PDEs in this model are written as follows:

$$\rho C\left(\frac{\partial }{\partial t}T+\mathbf{v}·\nabla T\right)-\nabla \left((4\sigma \delta T^3+K)\nabla T\right)+H(T-T_{\infty })=(T>T_p)QMA{e}^{-\frac{E_A}{RT}},$$

(7)

$$\frac{\partial }{\partial t}M=-(T>T_p)MA{e}^{-\frac{E_A}{RT}},$$

(8)

where $\sigma$ is the Stefan–Boltzmann constant, $\delta$ is the optical path length for radiation, and $T_p$ is the pyrolysis temperature.

A more rational change of variables is also incorporated for the non-dimensionalization of the equations, the Frank-Kamenetskii change of variables [25]. This change of variables involves setting reference values for temperature and fuel in which an equilibrium state can be assumed, as well as spatial and temporal variables, and permits elucidating the significant parameters. The reference values are the ambient temperature $T_{\infty }$ and the initial solid fuel load $M_0$. Then, the dimensionless temperature is now $u=(T-T_{\infty })/\epsilon T_{\infty }$, and the dimensionless solid fuel load is $c=M/M_0$. We kept again the notation for the spatial and temporal variables. The non-dimensional PDEs of this second model are

$$\frac{\partial }{\partial t}u+\mathbf{w}·\nabla u-\nabla (K\left(u\right)\nabla u)+\alpha u=f(u,c),$$

(9)

$$\frac{\partial }{\partial t}c=-\frac{\epsilon }{q}f(u,c),$$

(10)

where

$$\begin{array}{ccc}\hfill K\left(u\right)& =& \varkappa {\left(1+\epsilon u\right)}^3+1,\hfill \\ \hfill f\left(u,c\right)& =& (u>u_p)c{e}^{\frac{u}{1+\epsilon u}},\hfill \end{array}$$

$\alpha$ is the dimensionless natural convection coefficient, and it depends on the natural convection coefficient, the density, and the specific heat. The nonlinear function $K\left(u\right)$ in the diffusive term represents the local radiation. The parameter $\varkappa$ is the dimensionless inverse of the conductivity coefficient and depends on the Stefan–Boltzmann constant, the optical path length for radiation, the thermal conductivity, and the ambient temperature. $\epsilon$ is the dimensionless inverse of the activation energy and depends on the activation energy, the universal constant of gases, and the ambient temperature. The function $f(u,c)$ is given by an Arrhenius type law, where $u_p$ is the non-dimensional pyrolysis temperature, and the logic expression is equal to 1 if it is true and 0 if it is false. q is the non-dimensional reaction heat and depends on the heat of combustion, the initial fuel, and the specific heat. For a deep understanding of the dimensionless process, see [23].

Equations (9) and (10), together with appropriate initial and boundary conditions, provide challenging mathematical and numerical problems. Facing these types of mathematical questions is one of the challenges that should not be forgotten in the development of simulating models of processes as complex as wildfire spread. Some results about the existence and uniqueness of weak solutions of the non-convective version, that is the nonlinear reaction–diffusion problem, can be found in [23]. From the numerical point of view, the proposed scheme for the full version of the model is based on the use of the Mixed Finite-Element Method (MFEM), which admits discontinuities in the temperature, preserving the continuity of the flux through the inter-element boundaries. This allows the representation of high gradients in the solution corresponding to the fire front with strong temperature gaps. The scheme proposed uses the lowest-order Raviart–Thomas elements [26]. Finally, the convective term is solved by a splitting technique using Godunov’s method.

2.4. The Effect of Moisture Content: A Multivalued Operator in Enthalpy

The third improvement to the model was intended to reflect the effect of fuel moisture. The effect of the moisture content of the solid fuel was included through a multivalued maximal monotone operator relating enthalpy and temperature. The use of a multivalued operator was informed by the classical two-phase Stefan problem [27], and it was adapted to model the two well-defined phases in a wildfire combustion process [19]: the endothermic phase, which includes the dehydration of the solid fuel, and the exothermic phase, in which the flammable mixture from fuel pyrolysis begins to release energy. The Fuel Moisture Content (FMC) is one of the most-influential factors in fire spread, mainly through the process of heating and subsequently evaporating the water in the fuel, enabling it to attain combustion conditions. This process involves the consumption of the energy released by the adjacent combustion fuel and requires time, which reduces the fire Rate Of Spread (ROS) as moisture increases. The non-dimensional equations of this new model that appeared in [28] are, keeping the previous notation,

$$\frac{\partial }{\partial t}e+\mathbf{w}·\nabla e-\nabla (K\left(u\right)\nabla u)+\alpha u=f(u,c),$$

(11)

$$e\in G\left(u\right),$$

(12)

$$\frac{\partial }{\partial t}c=-\frac{\epsilon }{q}f(u,c),$$

(13)

where u and c are again the dimensionless temperature and mass fraction of solid fuel, respectively, and e is the dimensionless enthalpy.

The change of the Frank-Kamenetskii variable is maintained since an Arrhenius-type expression continues to appear in the reactive term. The diffusive term function $K\left(u\right)$, and the reactive term function $f(u,c)$ are the same as in the previous model. The multivalued operator is given by

$$G\left(u\right)=\left\{\begin{array}{ccc}u\hfill & \mathrm{if}& \hfill u<u_v,\\ [u_v,u_v+{\lambda }_v]\hfill & \mathrm{if}& \hfill u=u_v,\\ u+{\lambda }_v\hfill & \mathrm{if}& \hfill u_v<u<u_p,\\ [u_p+{\lambda }_v,u_p+{\lambda }_v+{\lambda }_p]\hfill & \mathrm{if}& \hfill u=u_p,\\ u+{\lambda }_v+{\lambda }_p\hfill & \mathrm{if}& \hfill u>u_p,\end{array}\right.$$

(14)

where $u_v$ is the dimensionless water evaporation temperature and $u_p$ is the dimensionless solid fuel pyrolysis temperature. The quantity ${\lambda }_v$ is the dimensionless evaporation heat related to the latent heat of evaporation ${\Lambda }_v$ and the fuel moisture content $M_v$ (kg of water/kg of dry fuel); ${\lambda }_p$ is the dimensionless pyrolysis heat related to the latent heat of pyrolysis ${\Lambda }_p$.

It is worth noting that, in the burnt area, the multivalued operator does not exactly represent the physical phenomena since water vapor is no longer in the porous medium. This inconvenience can be avoided by setting ${\lambda }_v=0$ and ${\lambda }_p=0$ in the burnt area.

The maximal monotone property of this multivalued operator allows the implementation of a numerical algorithm with good convergence properties based on the use of duality methods [29]. Given an exact perturbation of the multivalued operator, its properties, and an appropriate choice of the parameters, we can define the resolvent and its Yosida approximation, whereby the new nonlinear univalued operator equivalent to the multivalued one can be solved by a fixed-point iteration. For further details of how to numerically treat this multivalued operator, see [30,31], although Section 3.1 provides updated explanations about how the multivalued operator is numerically solved in the current model.

2.5. Non-Local Radiation: Some 3D Effects

In fact, the model in Section 2.4 was not implemented as-is, since another important improvement, the non-local radiation, was simultaneously incorporated. This is an essential improvement to reach the current version of the model, and it allows incorporating certain three-dimensional effects while maintaining the simplicity of a two-dimensional model.

In a first approximation, the idea of non-local radiation was introduced by means of a convolution operator, simultaneously maintaining the local radiation represented with the nonlinear term of the diffusive term of the energy equation; see [30]. However, based on the results, it was decided to simplify the energy equation and choose only non-local radiation [32]. This was a challenge from a computational point of view as it involved solving the radiation intensity equation in a 3D domain, the layer of air on the surface in which the fire develops.

The equations of the above models are defined over the surface S where the fire takes place, defined by the mapping:

$$\begin{array}{ccc}S:d\hfill & ⟼& {\mathbb{R}}^3\hfill \\ (x,y)\hfill & ⟼& (x,y,h(x,y\left)\right),\hfill \end{array}$$

(15)

where $h(x,y)$ is a known function representing the topography of the surface S.

In order to take into account some 3D effects, specifically non-local radiation, we shall consider the following 3D domain, representing the air layer of thickness $\delta$ on the surface S:

$$D=\left\{\right(x,y,z):x,y\in d,h(x,y)<z<h(x,y)+\delta \}.$$

The non-dimensional simplified equations defined on the surface S are

$$\frac{\partial }{\partial t}e+\mathbf{w}·\nabla e+\alpha u=r,$$

(16)

$$e\in G\left(u\right),$$

(17)

$$\frac{\partial }{\partial t}c=-g\left(u\right)c,$$

(18)

completed again with homogeneous Dirichlet boundary conditions, assuming that the surface S is large enough so that the fire does not reach the boundary during the simulation interval $t_{max}$, the initial conditions representing the initial fuel load and fire source.

The unknowns $e=\frac{E}{MC{T}_{\infty }}$, the dimensionless enthalpy, $u=\frac{T-T_{\infty }}{T_{\infty }}$, the dimensionless temperature of the solid fuel, and $c=\frac{M}{M_0}$, the mass fraction of solid fuel, are bidimensional variables defined in $S\times (0,t_{max})$. Note that we abandoned the Frank-Kamenetskii change of variable since an Arrhenius-type expression no longer appears in the energy equation. The physical quantities E, T, and M are the enthalpy, the temperature of solid fuel, and the fuel load, respectively; C is the heat capacity of solid fuel; $T_{\infty }$ is a reference temperature, namely the ambient temperature; $M_0$ is the initial fuel load.

The multivalued maximal monotone operator G is slightly different:

$$G\left(u\right)=\left\{\begin{array}{ccc}u\hfill & \mathrm{if}& \hfill u<u_v,\\ [u_v,u_v+{\lambda }_v]\hfill & \mathrm{if}& \hfill u=u_v,\\ u+{\lambda }_v\hfill & \mathrm{if}& \hfill u_v<u<u_p,\\ [u_p+{\lambda }_v,\infty ]\hfill & \mathrm{if}& \hfill u=u_p,\end{array}\right.$$

(19)

where, again, $u_v$, $u_p$, and ${\lambda }_v$ are the same as above. We simplified the model assuming that the maximum value that the dimensionless temperature of the solid fuel u reaches is the dimensionless solid fuel pyrolysis temperature $u_p$.

We must specify that, in this model, only the solid phase of the combustion process is considered: the mass fraction of solid fuel c is a dimensionless variable between 0 and 1, and as just mentioned, the maximum value of the dimensionless solid fuel temperature u is the dimensionless pyrolysis temperature $u_p$. Note that this is not the fire temperature, since the gaseous phase in this model is parameterized through the flame temperature $T_f$ and the flame height F in the radiation term, which at the moment are considered input data dependent on the type of fuel. In Section 2.6, we explain how to improve this parameterization.

We also assumed that the dimensionless solid fuel begins to be lost when it reaches this temperature, $u_p$. This means a simplification of the right-hand side of Equation (18), which represents the loss of solid fuel due to combustion. Now, $g\left(u\right)=0$ when $u<u_p$, and $g\left(u\right)$ is constant when $u=u_p$, where this constant is inversely proportional to the half-life time of the combustion of each type of fuel.

It remains to explain how the non-local radiation is computed, that is the r term on the right-hand side of Equation (16). Since radiation essentially comes from flames, we considered a simplified physical model in which the gases produced by pyrolysis burn above the fuel layer, producing a flame above that layer. The flame may eventually be tilted due to the wind or the slope of the ground and emits radiation, which reaches points ahead of the flame, heating the surrounding unburned fuel and allowing the fire to spread.

The term r describes the thermal radiation, which reaches the surface S from the flame above the layer

$$r=\frac{\left[t\right]}{MC{T}_{\infty }}R.$$

(20)

where $\left[t\right]$ is a time scale that appears during the dimensionless process. R represents the incident energy at a point $\mathbf{x}=(x,y,h(x,y\left)\right)$ of the surface S, due to radiation from the flame over the surface per unit time and per unit area, obtained by summing up the contribution of all directions $\Omega$, i.e.,

$$R\left(\mathbf{x}\right)={\int }_{\omega =0}^{2\pi }I(\mathbf{x},\Omega )\Omega ·\mathbf{N}d\omega ,$$

(21)

where $\omega$ is the solid angle and $\mathbf{N}$ is the unit normal vector to the surface S. We considered only the hemisphere above the fuel layer, and each contribution depends, among others, on the flame height F and the flame temperature $T_f$.

I is the total radiation intensity, described by the following differential equation:

$$\frac{dI}{ds}+a\left(s\right)I\left(s\right)=a\left(s\right)I_b\left(s\right),$$

(22)

where $I_b$ is the black body total radiation intensity, and it is governed by the Stefan–Boltzmann law:

$$I_b=\frac{\sigma }{\pi }{T}^4,$$

(23)

where $\sigma =5.6699\times {10}^{-8}\mathrm{w}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$ is the Stefan–Boltzmann constant and a is the radiation absorption coefficient inside the flame. The path s is within the 3D domain D reaching any point on the surface S and eventually passing through a flame. We assume that the temperature T reaches the flame temperature $T_f$.

Details on how to approximate Equation (22) and calculate $R\left(x\right)$ for different types of flame, vertical or tilted, can be found in [32]. We include here some aspects of the calculation of the radiation term for the case of a vertical flame in order to clarify some concepts that were introduced to reduce the computational time and to enhance the efficiency of the simulation process.

In our model, the radiation term R for a vertical rectangular flame has the form:

$$R\left(\mathbf{x}\right)=\frac{a\sigma T_f^4}{\pi }{\int }_{\Omega }\chi \left(\tilde{\mathbf{x}}\right)f(\mathbf{x}-\tilde{\mathbf{x}})d\tilde{A},$$

(24)

where $\chi$ is the characteristic function of the flames, i.e., $\chi \left(\mathbf{x}\right)=1$ if $T\left(\mathbf{x}\right)=T_p$ and $\chi \left(\mathbf{x}\right)=0$ otherwise, and $T_p$ represents the pyrolysis temperature.

We used the P1 finite-element method for the spatial discretization, so using the finite-element basis to represent the f function, the radiation term can be written for each node ${\mathbf{x}}_i$ of the mesh as follows:

$$R\left({\mathbf{x}}_i\right)=\frac{a\sigma T_f^4}{\pi }\sum _j{\mathrm{R}}_{ij}{A}_j\chi \left({\mathbf{x}}_j\right)$$

where $\mathrm{R}=(R_{ij}=f^j\left({\mathbf{x}}_i\right)),$ is the Radiation Matrix, representing the nodes of the finite element mesh, which are reached by the radiation emitted by each node, and it is computed only once (outside the time loop).

In practice, the radiation term is only calculated in a neighborhood of the fire front, which we call Active Nodes. At the beginning of the numerical process, we define a uniform and fine spatial mesh for spatial discretization, and for each time step $\Delta t$ of the time discretization, we define a set of Active Nodes formed by those nodes located inside a sufficiently large area close to the fire front (burning area); we solved the corresponding equations only in this set of nodes, as well as the calculation of the radiation term. Specifically, a node is Active if the dimensionless temperature in the node is positive, and the dimensionless solid fuel is $\ge 0.1$, or if the node belongs to the Radiation Molecule associated with a node fulfilling the previous conditions. The Radiation Molecule is a round set of nodes formed by the node itself and neighboring nodes that define the area affected by the radiation emitted by the node concerned. The Radiation Molecule was initially a set of 89 nodes formed by the node itself and 88 nodes surrounding it forming a circle of nodes of radius 5 nodes (see Figure 6 in [32]).

2.6. Flame Submodel

As discussed in the previous section, the calculation of non-local radiation strongly depends on the temperature and height of the flame, $T_f$ and F, respectively, which are initially considered model parameters depending on the type of fuel. These parameters have been improved by means of two submodels that delve into their physical meaning.

2.6.1. Flame Height Submodel

Flame angle and flame length are determining factors in the radiative heat transfer of a flame. Flame geometry plays a key role in determining the rate of spread of fires where radiation is the dominant heat transfer mechanism to the unburned fuel [33]. The joint effects of wind speed and slope percentage affect the flame geometry. We focused on the flame height F, understanding that this is the height of the flame above the vertical, relating the flame length $L_f$ and flame angle ${\theta }_f$, $F=L_{f}sin{\theta }_f$. The flame height, and not the flame length, is one of the input parameters determining the non-local radiation calculation in our model. Data from [34] and the sensitivity analysis of the model performed in [35] led us to propose the following flame height sub-model:

$$F=(F_H+F_v{\left|v\right|}^{1/2})(1+F_s{s}^2)$$

(25)

where $F_H$ is a flame-height-independent parameter, $F_v$ is a wind correction factor, $F_s$ is a slope correction factor, $\left|v\right|$ is the wind strength, and s represents the slope at each point on the surface. The first factor in Equation (25) represents the correction of the flame height due to the wind. This expression is derived from the observation of the experimental curves for different fuels gathered in [34], where the increase in the height of the flame due to the wind responds to such a function. $F_H$ stands for zero wind, and the correction coefficient $F_v$ was added to experimentally adjust the different behaviors for each type of fuel by a least-squares method. The second factor in Equation (25) represents the correction of the flame height according to surface slope. When there is no slope or wind, the flame is vertical, and its height is equal to its length; as deduced from [36], a relationship can be derived between the flame height and the vertical gas velocity inside it. However, when there is a slope, the following proportionality factor must be considered:

$$\frac{1}{{cos}^2\alpha }=1+{\left(\frac{\partial h}{\partial x}\right)}^2+{\left(\frac{\partial h}{\partial y}\right)}^2=1+s^2$$

where $\alpha$ is the angle between the horizontal plane and tangential plane to the surface S on each point and $h=h(x,y)$ is the surface height, so s is the slope at each point on the surface S. Again, a correction factor $F_s$ was added to adjust data from [34] in the sense of least-squares.

2.6.2. Flame Temperature Submodel

In our simplified fire spread model, the gaseous fuel is parameterized by the flame temperature in the non-local radiation calculation. We can improve this input parameter by writing the generic Equation (1) for the temperature inside the flame $T_f$ and assuming some simplifications. We assumed that this temperature does not vary inside a stabilized flame. We also considered that heat losses inside the flame are mainly due to local radiation, representing this term as $h_r(T_f^4-T_{\infty }^4)$, and that the heat source due to gas combustion is represented by $QM$, where M represents the solid fuel and Q is the heat generated or absorbed per unit mass of solid fuel. Since we wanted to represent the heat generated, that is when the process is exothermic, we added the term $(t>t_p)$, which is 0 when the process is endothermic and 1, when it is exothermic. Finally, the temperature equation inside a stabilized flame can be written as follows:

$$h_r(T_f^4-T_{\infty }^4)=(t>t_p)QM$$

(26)

Rearranging the previous expression, we obtained the following value for the flame temperature:

$$T_f={(T_{\infty }^4+(t>t_p)\frac{Q}{h_r}M)}^{1/4}$$

(27)

If an estimate of the maximum flame temperature $T_{f,max}$ is available, it can be used as the input in Equation (26) to obtain an approximation of $Q/h_r$ so that the following approximation of the flame temperature in terms of solid fuel is finally obtained:

$$T_f={(T_{\infty }^4+(t>t_p)(T_{f,max}^4-T_{\infty }^4)M)}^{1/4}$$

(28)

2.7. Fire-Spotting

The spot fire ignition of a wildland fire by hot sparks and firebrands is an important fire ignition pathway by which wildfires are started and may propagate. Once the wildfire has been ignited and grows, it can spread rapidly through ember spotting, where pieces of burning material are lofted by the plume of the fire and then transported forward by the wind, landing where they can start spot fires downwind beyond the main fire’s direct ignition zone. These secondary fires frequently cause hazardous situations for firefighters and contribute to the increase in the rate of fire spread. There are decisive factors affecting the generation of fire-spotting and the firebrands’ landing patterns and flight paths. Some are macro-scale factors such as the atmospheric stability [37] and others are meso-scale factors such as the flame geometry [38]. Nevertheless, fire-spotting has a strong random component, so available fire-spotting models are inherently probabilistic.

In [18], we proposed a preliminary fire-spotting sub-model that follows the characteristics of our model based on the ideas of the RandomFront 2.3 in [39], as a random heat contribution added to the right-hand side of the energy conservation equation. The heat contribution due to fire-spotting is written in terms of the distance of firebrand distribution $\varphi (\ell )$:

$$q=F_q\times N_q\times {\Delta }_t\times \varphi (\ell )$$

(29)

where $F_q$ is a factor transforming the probability density function $\varphi (\ell )$ into energy, $\Delta t$ is the time step of the time discretization used to numerically solve the PhyFire model equations, and $N_q$ is a fire-spotting index introduced to significantly reduce the computational effort of the fire-spotting module. In practice, this computation is performed every $N_q\times \Delta t$, instead of every $\Delta t$.

The firebrand landing distance ℓ can be assumed to follow a log-normal distribution [40]:

$$\varphi (\ell )=\frac{1}{\ell \sigma \sqrt{2\pi }}exp(-\frac{{(ln(\ell )-\mu )}^2}{2{\sigma }^2})$$

(30)

where $\mu$ and $\sigma$ are the mean and standard deviation of the logarithm of the variable ℓ, that is $ln\ell \sim N(\mu ,\sigma )$.

The set of possible firebrand emitter nodes is computed as a subset of the fire front. Next, the set of possible receiving nodes of the firebrand from each emitting node is computed as a random subset such that the main direction of firebrand propagation is the wind direction, and the firebrands’ landing distance depends on the mean distance of firebrand landing. The details of these calculations can be found in [18].

3. Current Model and Numerical Algorithm

All the improvements and sub-models described in the previous sections have been included in the current fire propagation model, PhyFire, as well as certain computational improvements, which we describe next.

3.1. Model Description

The surface S where the fire takes place is defined by the mapping described by the expression (15). The size of this surface S and the study time interval $(0,t_{max})$ were selected so that the fire does not reach the boundary in time $t_{max}$, to assume Dirichlet homogeneous boundary conditions to complete the set of PDEs that governs the PhyFire model. The latest version of these non-dimensional PDEs are as follows:

$$\begin{array}{c}\hfill {\partial }_{t}e+\beta \mathbf{v}·\nabla e+\alpha u=r+q\mathrm{in}S\times (0,t_{max}),\end{array}$$

(31)

$$\begin{array}{c}\hfill e\in G\left(u\right)\mathrm{in}S\times (0,t_{max}),\end{array}$$

(32)

$$\begin{array}{c}\hfill {\partial }_{t}c=-g\left(u\right)c\mathrm{in}S\times (0,t_{max}),\end{array}$$

(33)

where the unknowns are the following variables defined in $S\times (0,t_{max})$: the dimensionless enthalpy e, the dimensionless solid fuel temperature u, and the solid fuel mass fraction c. The relationship between these dimensionless variables and the corresponding physical quantities is given in

Table 1. Dimensionless unknowns and related physical quantities.

Physical Variable	Symbol	Units	Dimensionless Variable
Enthalpy	E	J m${}^{-2}$	$e=E/MC{T}_{\infty }$
Solid fuel temperature	T	K	$u=(T-T_{\infty })/T_{\infty }$
Solid fuel load	M	kg m${}^{-2}$	$c=M/M_0$

.

The physical quantities that appear in Table 1 are the enthalpy E, the solid fuel temperature T, and the fuel load M. Their relationship depends on the reference temperature $T_{\infty }$, the heat capacity of the solid fuel C, and the maximum solid fuel load $M_0$. The reference temperature $T_{\infty }$ is related to the ambient temperature, measured far enough from the fire front to ensure that $T\ge T_{\infty }$, whereby $u\ge 0$. The heat capacity of the solid fuel C and maximum solid fuel load $M_0$ both depend on the fuel type.

Table 2. Fuel-type-dependent input variables.

Input Variable	Symbol	Units
Heat capacity	C	J K${}^{-1}$ kg${}^{-1}$
Maximum initial fuel load	$M_0$	kg m${}^{-2}$
Fuel moisture content	$M_v$	kg water/kg fuel
Maximum flame temperature	$T_{f,max}$	K
Pyrolysis temperature	$T_p$	K
Combustion half-life	$t_{1/2}$	s
Flame length independent factor	$F_H$	m
Flame length wind correction factor	$F_v$	m${}^{1/2}$s${}^{1/2}$
Flame length slope correction factor	$F_s$	−

lists all the variables that depend on the fuel type.

The PhyFire model also depends on the three model parameters listed in

Table 3. Model parameters.

Parameter	Symbol	Units
Natural convection coefficient	H	J s${}^{-1}$m${}^{-2}$K${}^{-1}$
Convective term factor	$\beta$	−
Mean absorption coefficient	a	m${}^{-1}$

. These three parameters are related to the heat transfer terms in Equation (31): H in the natural convection term $\alpha u$, $\beta$ in the convective term $\beta \mathbf{v}·\nabla e$, and a in the radiation term r.

We completed the problem with the following initial conditions:

$$\begin{array}{c}\hfill u(x,y,0)=u_0(x,y)\mathrm{in}S,\end{array}$$

(34)

$$\begin{array}{c}\hfill c(x,y,0)=c_0(x,y)\mathrm{in}S.\end{array}$$

(35)

which represent, respectively, the focus of the fire, or eventually an intermediate perimeter if it is decided to restart the simulation from an instant later than the initial one, and the initial distribution of the fuel, including any possible firebreaks.

We should point out that the data for topography (function $h(x,y)$ defining S), solid fuel distribution (initial value of M), and fuel type distribution (to spatially determine the fuel-type-dependent input variables) were obtained from the cartography generated particularly for each case. This process was automated by developing the necessary cartographic database [17] and the online GIS interface [18]. Section 4 details some aspects related to the cartographic and meteorological information necessary to apply the PhyFire model in real cases.

The zero-order term $\alpha u$ in Equation (31) represents the energy lost by natural convection in the vertical direction, as in the initial models. Here, parameter $\alpha$ depends on the natural convection coefficient H through the following expression:

$$\alpha =\frac{H\left[t\right]}{MC}$$

(36)

where $\left[t\right]$ is a time scale.

The convective term, $\beta \mathbf{v}·\nabla e$, represents the energy convected by the gas pyrolyzed through the elementary control volume. $\mathbf{v}$ is the wind velocity re-scaled by a correction factor $\beta$. The idea behind this parameter is how the one-phase model, PhyFire (mainly a solid phase model), can be understood as a simplification of a gas–solid two-phase model where we retained the enthalpy transported in the gas phase. The enthalpy transported in the gas phase is retained, and $\beta$ represents the fraction of transported heat that must be taken into account in Equation (31). For typical values of the heat capacity for the air and for the fuel (wood), $\beta$ can be estimated to have an order of magnitude of ${10}^{-2}$ inside the flame and ${10}^{-4}$ away from the flame. See [41] for more details.

Furthermore, we have to mention that the surface wind velocity $\mathbf{v}$ can be considered as the given data or can be computed, for example, by means of the high-definition wind field model HDWind specifically developed to be used coupled with the PhyFire model (see [14,15]).

As we explained in Section 2.4, the multivalued Equation (32) models the influence of moisture content and depends on the fuel moisture content $M_v$ and solid fuel pyrolysis temperature $T_p$. The dimensionless enthalpy e is, thus, an element of the multivalued maximal monotone operator G defined by Equation (19), where, remember, $u_v$ is the non-dimensional water evaporation temperature and $u_p$ is the non-dimensional solid fuel pyrolysis temperature. It remains to specify the expression defining the non-dimensional evaporation heat ${\lambda }_v$:

$${\lambda }_v=\frac{M_v{\Lambda }_v}{C{T}_{\infty }}$$

(37)

where the latent heat of water evaporation ${\Lambda }_v=2.25\times {10}^6$$Jk{g}^{-1}$ and the fuel moisture content $M_v$ (kg of water/kg of dry fuel) is one of the fuel-type-dependent input variables in Table 2. We should clarify that the fuel moisture content depends not only on the fuel-type, but also on other factors that vary both spatially and temporally. Fuel moisture content is driven by several weather factors, such as long-term droughts, recent rains, or immediate temperature and relative humidity [42]. Topographic aspect, i.e., the direction that the land faces, is another factor that affects the fuel moisture content [43]. Depending on the landscape position, there are areas that receive more radiation, will dry more quickly, and will, therefore, have lower fuel moisture. One of the latest enhancements to the PhyFire model is the ability to use space-varying fuel moisture content data, replacing the fuel-type-dependent input variable $M_v$ with a raster containing the spatial distribution of the $M_v$ values, using remote sensing data.

The non-local radiation term r in the right-hand side of Equation (31) represents the radiation from the flames above the surface where the fire takes place, as we explained in Section 2.5. This term depends on the incident energy due to radiation R received in each point of the surface S given by Equation (20). The incident energy R should be computed by adding the contribution of all directions on the hemisphere above the fuel layer following the Equation (21), where the total radiation intensity I is given by the differential Equation (22). Each incident energy contribution depends on the flame temperature $T_f$, which is approximated by Equation (28) once dimensionless, and the flame height F, which is calculated depending on the wind, slope, and fuel type, according to Equation (25).

The term q on the right-hand side of of Equation (31) stands for the random heat contribution due to fire-spotting, explained in Section 2.7.

Equation (33) represents the loss of solid fuel due to combustion, where

$$g\left(u\right)=\left\{\begin{array}{ccc}0\hfill & if& \hfill u<u_p,\\ \gamma \hfill & if& \hfill u=u_p,\end{array}\right.$$

(38)

i.e., there is no loss of solid fuel if the temperature is below the pyrolysis temperature, and it remains constant when the temperature of pyrolysis is reached. This constant value is inversely proportional to the solid fuel half-life time $t_{1/2}$, of the combustion of each type of fuel:

$$\gamma =\frac{ln2\left[t\right]}{t_{1/2}}.$$

(39)

3.2. Numerical Method

The PhyFire code used here incorporates a new numerical scheme based on the use of P1 finite-element approximation on a regular mesh for spatial discretization, as in the previous schemes. However, a new predictor–corrector finite difference scheme was developed for time discretization. The predictor step is an Euler semi-implicit scheme, and the corrector step is a modified Crank–Nicolson scheme. The computational cost was reduced by defining the set of Active Nodes, as we explain at the end of Section 2.5, solving the equations only in the neighborhood of the fire front and adapting the numerical scheme and the corresponding code to parallel computing [32].

Given the initial values $u^0$ and $c^0$ defined by the initial conditions, we set the value of the initial enthalpy for each node i of the spatial discretization depending on the initial nondimensional temperature $u^0$ as follows:

$$e_i^0=\left\{\begin{array}{ccc}{u}_i^0\hfill & \mathrm{if}& \hfill u_i^0\le u_v,\\ u_i^0+{\lambda }_v\hfill & \mathrm{if}& \hfill u_i^0>u_v,\end{array}\right.$$

(40)

where, as we have mentioned before, $u_v$ is the non-dimensional water evaporation temperature and ${\lambda }_v$ is the dimensionless evaporation heat. As we explained in the description of the multivalued operator G, we assumed that the maximum value that the solid fuel temperature reaches is the pyrolysis temperature, because the gaseous phase is parameterized, so the only phase change in the temperature of the solid fuel corresponds to the water evaporation.

Given the values of the unknowns $u^n$, $c^n$, and $e^n$ at time step n, we computed $u^{n+1}$, $c^{n+1}$, and $e^{n+1}$ by means of the following steps:

• Build the set of Active Nodes.
• Compute the Radiation Heat.
• Prediction step: semi-implicit Euler method.
• Update the set of Active Nodes.
• Update the Radiation Heat.
• Correction step: modified Crank–Nicolson method.

Lets us remember that a node i is considered an Active Node if $u_i>0$ and $c_i\ge 0.1$ or if it belongs to the Radiation Molecule associated with a node fulfilling the previous condition. The number of nodes defining the Radiation Molecule was modified in order to keep its real size independent of the spatial discretization size. Specifically, a number of nodes of the radius of the round radiation molecule was fixed for each mesh size. The mesh size of the spatial discretization varied depending on the precision level, which currently varied from precision level 0 corresponding to 50 m cell size, to simulate large fires, to Precision Level 5, corresponding to a 2.5 m cell size, to simulate small fires. We also contemplated much thinner meshes, on the order of centimeters, to simulate laboratory experiments. In any case, the input raster files were adapted to the FEM precision level.

3.2.1. Convection Step

We must begin by explaining the total discretization of the convective term of (31), which was carried out in each step of the predictor–corrector scheme:

$${\partial }_{t}e+\beta \mathbf{v}.\nabla e\approx \frac{1}{\Delta t}(e^{n+1}-{\overline{e}}^n)$$

where ${\overline{e}}^n=e^n\circ X^n$ and $X^n\left(\mathbf{x}\right)=X(\mathbf{x},t^{n+1},t^n)\approx \mathbf{x}-\beta \mathbf{v}\Delta t$ is the position at time $t^n$ of the particle that is at position $\mathbf{x}$ at time $t^{n+1}$ [44].

3.2.2. Predictor Step

The discrete equations in the burning area in the predictor step respond to a semi-implicit Euler scheme:

$$\begin{array}{cc}\hfill & \frac{e^{n+1/2}-{\overline{e}}^n}{\Delta t}+\alpha u^{n+1/2}=r^n+q^n,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & e^{n+1/2}\in G\left(u^{n+1/2}\right),\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & \frac{c^{n+1/2}-c^n}{\Delta t}=-g\left(u^{n+1/2}\right)c^{n+1}.\hfill \end{array}$$

(43)

where $e^{n+1/2}$, $u^{n+1/2}$, and $c^{n+1/2}$ mean the predicted value of $e^{n+1}$, $u^{n+1}$, and $c^{n+1}$, respectively. The basic idea is to treat the linear term implicitly. The heat contribution of non-local radiation r and fire-spotting q depends on the dimensionless temperature u and solid fuel mass c, so the non-local radiation term $r^n$ and fire-spotting term $q^n$ are evaluated explicitly at time $t^n$. Even so, Equations (41)–(43) continue to be nonlinear due to the multivalued operator G. However, the solution to this problem can be reduced to explicit calculations.

The multivalued operator in Equation (42) is maximal monotone, and hence, its resolvent $J_{\mu }={(Id+\mu G)}^{-1}$ is a well-defined univalued operator for any $\mu >0$. Moreover, the Yosida approximation [29] of G, $G_{\mu }=\frac{Id-J_{\mu }}{\mu }$ is a Lipschitz operator, and the inclusion in Equation (42) is equivalent for all $\mu >0$ to the following equation:

$$e^{n+1/2}=G_{\mu }(u^{n+1/2}+\mu e^{n+1/2}),$$

(44)

or

$$u^{n+1/2}=J_{\mu }(u^{n+1/2}+\mu e^{n+1/2}).$$

(45)

Now, rearranging Equation (41):

$$u^{n+1/2}+\frac{1}{\alpha \Delta t}{e}^{n+1/2}=\frac{1}{\alpha \Delta t}{\overline{e}}^n+\frac{1}{\alpha }(r^n+q^n).$$

(46)

and taking $\mu =1/(\alpha \Delta \tau )$ in Equation (45), we obtain

$$u^{n+1/2}=J_{1/\alpha \Delta t}\left(\frac{1}{\alpha \Delta t}{\overline{e}}^n+\frac{1}{\alpha }(r^n+q^n)\right)$$

(47)

which it is equivalent to solving

$$(\alpha \Delta tId+G)u^{n+1/2}\ni {\overline{e}}^n+\Delta t(r^n+q^n).$$

(48)

Thus, denoting ${\ae }^n={\overline{e}}^n+\Delta t(r^n+q^n)$, the value of $u^{n+1/2}$ is given by

$$\begin{array}{ccc}\hfill \frac{{\ae }^n}{1+\alpha \Delta t}& \mathrm{if}& {\ae }^n<(1+\alpha \Delta t)u_v,\hfill \\ \hfill u_v& \mathrm{if}& {\ae }^n\in [(1+\alpha \Delta t)u_v,(1+\alpha \Delta t)u_v+{\lambda }_v],\hfill \\ \hfill \frac{{\ae }^n-{\lambda }_v}{1+\alpha \Delta t}& \mathrm{if}& {\ae }^n\in [(1+\alpha \Delta t)u_v+{\lambda }_v,(1+\alpha \Delta t)u_p+{\lambda }_v],\hfill \\ \hfill u_p& \mathrm{if}& {\ae }^n\in [(1+\alpha \Delta t)u_p+{\lambda }_v,\infty ].\hfill \end{array}$$

(49)

Once $u^{n+1/2}$ has been obtained, we calculate $e^{n+1/2}$ and $c^{n+1/2}$ explicitly from Equation (41) and (43), respectively,

$$\begin{array}{c}\hfill e^{n+1/2}={\overline{e}}^n-\alpha \Delta t{u}^{n+1/2}+\Delta t(r^n+q^n),\end{array}$$

(50)

$$\begin{array}{c}\hfill c^{n+1/2}=\frac{c^n}{1+\Delta tg\left(u^{n+1/2}\right)}.\end{array}$$

(51)

Notice that Equations (47), (50), and (51) can be solved simultaneously in all active nodes, and hence, parallel computation can be used to reduce the computational cost. In actual fact, the loop over all actives nodes to compute $u^{n+1/2}$, $e^{n+1/2}$, and $c^{n+1/2}$ was parallelized using the API OpenMP [45].

In the burned area, only Equation (41) needs to be solved. Since there is no fuel in this area and there are no changes in enthalpy or fire-spotting, Equation (41) simplifies to

$$\begin{array}{c}\hfill \frac{u^{n+1/2}-u^n}{\Delta t}+\alpha u^{n+1/2}=r^n\end{array}$$

(52)

which can be solved explicitly:

$$\begin{array}{c}\hfill u^{n+1/2}=\frac{u^n+\Delta t{r}^n}{1+\alpha \Delta t}\end{array}$$

(53)

3.2.3. Corrector Step

The discrete equations in the burning area in the corrector step respond to a Crank–Nicolson scheme:

$$\begin{array}{cc}\hfill & \frac{e^{n+1}-{\overline{e}}^n}{\Delta t}+\alpha \frac{u^{n+1}+u^n}{2}=\frac{r^{n+1/2}+r^n}{2}+\frac{q^{n+1/2}+q^n}{2}\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill & e^{n+1}\in G\left(u^{n+1}\right),\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill & \frac{c^{n+1}-c^n}{\Delta \tau }=-g\left(\frac{u^{n+1}+u^n}{2}\right)\frac{c^{n+1}+c^n}{2}\hfill \end{array}$$

(56)

where we approximated $r^{n+1}$ and $q^{n+1}$ by $r^{n+1/2}$ and $q^{n+1/2}$, respectively, computed in terms of $u^{n+1/2}$ and $c^{n+1/2}$, the estimations obtained in the prediction steps, as at this point, $u^{n+1}$ and $c^{n+1}$ are not known. As in the prediction step, rearranging Equation (54), we have:

$$u^{n+1}+\frac{2}{\alpha \Delta t}{e}^{n+1}=\frac{2}{\alpha \Delta t}{\overline{e}}^n-u^n+\frac{r^{n+1/2}+r^n}{\alpha }+\frac{q^{n+1/2}+q^n}{\alpha }.$$

(57)

As in the predictor step, taking $\mu =\frac{2}{\alpha \Delta t}$, the multivalued Equation (55) can be written as

$$\begin{array}{c}\hfill u^{n+1}=J_{\frac{2}{\alpha \Delta t}}\left(\frac{2}{\alpha \Delta t}{\overline{e}}^n-u^n+\frac{r^{n+1/2}+r^n}{\alpha }+\frac{q^{n+1/2}+q^n}{\alpha }\right)\end{array}$$

(58)

which it is equivalent to solving

$$(\frac{2}{\alpha \Delta \tau }Id+G)u^{n+1}\ni \left(\frac{2}{\alpha \Delta t}{\overline{e}}^n-u^n+\frac{r^{n+1/2}+r^n}{\alpha }+\frac{q^{n+1/2}+q^n}{\alpha }\right).$$

(59)

Denoting now ${\ae }^n=\frac{2}{\alpha \Delta t}{\overline{e}}^n-u^n+\frac{r^{n+1/2}+r^n}{\alpha }+\frac{q^{n+1/2}+q^n}{\alpha }$, the value of $u^{n+1}$ is given by

$$\begin{array}{ccc}\hfill \frac{{\ae }^n}{1+\frac{\alpha \Delta t}{2}}& \mathrm{if}& {\ae }^n<(1+\frac{\alpha \Delta t}{2})u_v,\hfill \\ \hfill u_v& \mathrm{if}& {\ae }^n\in [(1+\frac{\alpha \Delta t}{2})u_v,(1+\frac{\alpha \Delta t}{2})u_v+{\lambda }_v],\hfill \\ \hfill \frac{{\ae }^n-{\lambda }_v}{1+\frac{\alpha \Delta t}{2}}& \mathrm{if}& {\ae }^n\in [(1+\frac{\alpha \Delta t}{2})u_v+{\lambda }_v,(1+\frac{\alpha \Delta t}{2})u_p+{\lambda }_v],\hfill \\ \hfill u_p& \mathrm{if}& {\ae }^n\in [(1+\frac{\alpha \Delta t}{2})u_p+{\lambda }_v,\infty ].\hfill \end{array}$$

(60)

Again, once $u^{n+1}$ has been obtained, we can update the enthalpy $e^{n+1}$ from Equation (54) and the fuel $c^{n+1}$ from Equation (56):

$$\begin{array}{c}\hfill e^{n+1}=e^n-\alpha \Delta t\frac{u^{n+1}+u^n}{2}+\Delta t\frac{r^{n+1/2}+r^n}{2}+\frac{q^{n+1/2}+q^n}{2},\end{array}$$

(61)

$$\begin{array}{c}\hfill c^{n+1}=\frac{1-\frac{\Delta t}{2}g\left(\frac{u^{n+1}+u^n}{2}\right)}{1+\frac{\Delta t}{2}g\left(\frac{u^{n+1}+u^n}{2}\right)}{c}^n.\end{array}$$

(62)

Once again, Equations (58), (61), and (62) can be solved simultaneously in all active nodes, so again, the parallel calculation allowed us to reduce the computational cost in the correction step.

In the burned zone, only (54) needs to be taken into account, which can be simplified as well and solved explicitly:

$$\begin{array}{c}\hfill u^{n+1}=\frac{(1-\frac{\alpha \Delta t}{2})u^n+\frac{\Delta \tau }{2}(r^{n+1/2}+r^n+q^{n+1/2}+q^n)}{1+\frac{\alpha \Delta t}{2}}\end{array}$$

The computational cost was reduced by defining the Active Nodes, solving the equations only for the fire front, computing the Radiation Matrix once out of the time loop, differentiating the burning zone from the burned zone, and adapting the numerical scheme and the corresponding code to parallel computing. We also introduced the fire-spotting index $N_q$ to significantly reduce the computational cost of the set of possible firebrand receiver nodes, so that the term $q^n$ is not calculated in all the time steps, but every $N_q$ steps, although for greater clarity, we detailed the calculation of this term in each time step. Improving computational efficiency is crucial to obtain simulation results in times that are significantly lower than the real-time evolution of a fire. However, it is also of great importance to address other objectives in the design and improvement of this type of complex models, such as the sensitivity analysis or the parameter adjustment, which require a large number of model evaluations, as well as to improve simulation results through data assimilation. All these aspects have been analyzed for the PhyFire model in different previous works [32,35,41].

The PhyFire model was developed under C++, using Neptuno++, a finite-element toolbox mainly developed by L. Ferragut [46], and has the advantages of the shared memory parallel paradigm using the OpenMP API.

4. GIS Integration

The ultimate goal of this kind of environmental model was to provide the potential end-user with a simple, intuitive, and easy-to-use tool. The key to achieve this was to automate the pre-processing of the input data and the post-processing of the simulation results, showing them in a user-friendly interface. Environmental event simulation requires coupling geospatial data with the simulation models, and this approach calls for a geospatial data management system to handle the diverse resources of information needed as the input data, as well as to display the simulation results. The Geographical Information System (GIS) is one of the geospatial information technologies that can be used to store, manage, and analyze geographical data about natural factors contributing to wildfires, and therefore, it is a suitable tool to integrate a wildfire simulation model.

Different attempts have been made to integrate the PhyFire model in a GIS system: as an add-in for the desktop software ArcMap and as a GIS solution based on a web server. In [17], several scripts were developed to automate geospatial data acquisition and pre-processing, model execution, and visualization of the simulation results on a base map. This implementation used the Python language and Esri’s ArcPy library to extend the functionality of ArcMap 10.4. Although fully functional, integrating the PhyFire model into ArcGIS makes it difficult for potential usersto obtain the model, as it requires them to have the appropriate ArcGIS Desktop software license. To overcome these disadvantages, make the PhyFire model a widely accessible tool, and reach a broader audience, a web platform was developed in [18]. The platform http://sinumcc.usal.es (accessed on 30 January 2023) allows uploading fire simulation data and the pre-processing, processing, and visualization of the simulation results throughout the Spanish territory. The web platform involves two large modules: the Sinumcc Web Platform (frontend) itself and the Sinumcc GIS Server (backend). The Sinumcc Web Platform is used to collect data intuitively, semi-automatically, and visually, as well as to display the results of the simulations. The Sinumcc GIS Server is the module that processes the spatial data for the area selected by the users and performs the simulations based on the data collected from the Sinumcc Web Platform. The web platform was developed using current communication and data processing technologies, such as API REST, JSON, and ArcGIS Server. Both solutions need a base map to help the user identify the different points of reference of the study domain: both the Spanish Topographic Base Map provided by the Instituto Geográfico Nacional de España (IGN) via a Web Map Service [47] and the Topographic base map from Esri were used for this purpose [48]. Currently, we are working on a profound improvement of this web platform to make it more accessible and expandable through the use of free software. Maintaining the philosophy of a frontend and backend (client–server) connected through web services based on a REST API, the web platform is being redesigned through a module-based design. Technologically, the new platform will be based on the use of a No-SQL database for user management, a spatial data model of global scope based on the World Geodetic System WGS84 [49] geographic coordinate system, enabling better access to geospatial information, the free software library GDAL/OGR [50], eliminating its dependence on the ArcGIS Server, the framework Angular 2+ [51], and OpenLayers as the map viewer.

For any of the GIS solutions adopted, it is essential to gather all the necessary geospatial information, as well as meteorological data, to feed the PhyFire model. For this purpose, a geodatabase was developed for the Spanish territory, which contains three maps with information adapted to the requirements of the PhyFire model.

The first map integrated into our geodatabase is the topographic map, which contains the height h at each point of the surface S, defined by the mapping given by Equation (15). This information can be found in any Digital Elevation Model (DEM), and we initially selected the DEM also published by the IGN via a Web Coverage Service (WCS) [52] with a resolution of 25 m, currently deprecated. Today, cartographic products with the same purpose and resolutions of 2 [53] and 5 [54] meters are available. The topographic map with a 5 m resolution is the most-appropriate option for PhyFire’s features.

The second map collects all the elements involving the function of either artificial or natural fuelbreaks that affect the fire spread. These were initially extracted from the Spanish Land Cover Information System (SIOSE) [55] by selecting all areas where a fire cannot occur (barren land, water bodies, transport infrastructures, etc.), providing zero-load fuel data for the model. It should be noted that forest fires simulated with the PhyFire model can cross the area associated with these elements thanks to the modeling of the radiation process. Today, new versions of the SIOSE product [56] together with other highly detailed products can be used for the same purpose, such as the National Topographic Database at a 1:25,000 scale (BTN25) [57], or the collaborative project OpenStreetMap [58]. It is possible to include temporal scenarios, especially those related to the different states of the crops, by using satellite data and the estimation of different spectral indexes such as the NDVI [59] or the NDBR [60].

The third map gathers all the information related to the fuel type, which is especially important for the PhyFire model, as most of its input variables depend on the fuel type (see Table 2). The PhyFire model can be tuned to work with different fuel classification systems, such as the Northern Forest Fire Laboratory (NFFL) system [61], the Fire Behavior Fuel Models (FBFM) [62], and of course, the Mediterranean-European Prometheus system [63], or even monitored data whenever available. In Spain, depending on the region and/or year in which the fire to be simulated occurred, the fuel type map used by PhyFire was processed by assembling the Spanish Forestry Map 1:25,000 (MFE25) [64] and the Spanish Forestry Map 1:50,000 (MFE50) [65]. Some limitations of this type of information are its lack of update, the lack of information in some areas, the non-consideration of the seasonal nature of forest fuel, or the subjective interpretation in the field work when this information is gathered. It is possible to solve these problems by using satellite data [63] for the specific instants when the simulation is going to be carried out and algorithms of machine learning [66] to keep the fuel classification updated.

The pre-processing of the geospatial input data and the post-processing of the simulation results can slow down the simulation process due to the high amount of information to be handled. New improvements have been developed to further enhance the efficiency of the simulation system embedding part of the pre- and post-processing functionalities through the integration of GDAL/OGR [50] utilities in the original code of the PhyFire model.

5. Real Case

In this section, we present the simulation of a real wildfire with the current PhyFire version, which includes the flame temperature and height submodels, as well as fire-spotting. We selected this wildfire for two reasons: firstly because we have all the relevant information of this fire, which was documented by the coordinator of the fire-suppression operations in [67], and, secondly, because the extensive photo-guide [34] provides us suitable information about the characteristics of the fuel types of the simulation area.

This wildfire occurred in Osoño, Ourense province, in the autonomous region of Galicia in north-western Spain, on 17 August 2009. The fire ignited at 3:45 p.m. local time near the motorway bridge. The firefighting team had failed to stabilize the fire by 11:00 p.m. on the same day, but brought it under control at 3:45 a.m. on the following day, finally extinguishing it at 9:10 p.m. on 18 August. The fire burned 224 ha:185 ha of forest area (83 ha were tree-covered interspersed with heath) and 39 ha of agricultural area.

The fuel data were obtained from the IFN4 with the fuel type distribution according to the NFFL classification [61], adapted to Spanish forestry by the Nature Conservation Institute [68]. The initial burnt area was covered mainly with Pinus pinaster corresponding to Model 7 (inflammable brush); the middle area was covered with Model 6 (dormant brush) and, to a lesser extent, Model 2 (timber grass). The end burnt area was covered with diverse fuel types, mainly Model 5 (brush).

Figure 2 shows the fuel type distribution, as well as the existing fuelbreaks, roads, rivers, and particularly interesting, a wide fuelbreak (20 m wide) along the crests that the actual fire crossed, as did the simulation.

As can be seen in the contour lines of Figure 2, for the first hours, the fire spread over an uneven surface, with positive and negative discontinuous slopes, with watersheds and river basins; for the final part, although the altitude is higher, the surface is relatively flat, and the terrain is more accessible for firefighter teams.

The firefighting team carried out a series of actions to prevent the fire from reaching the municipality of Osoño, establishing a set of firebreaks by widening the existing roads on the southern flank of the fire. Figure 2 also shows the three firebreaks (blue lines) the firefighters made. We did not consider in the simulation the actions of the fire truck water tanks on both flanks during the last part of the fire.

We simulated the first and most-complicated hours of this real fire, specifically 4 h and 15 min, updating the wind data (wind speed and direction) every hour and providing graphic results of the burnt area and fire front every 15 min. The simulation area was a rectangle measuring 3320 m × 2745 m; the cartographic spatial resolution was 5 m; and the finite-element mesh resolution was 7.5 m. The whole simulation involved about 150 s of computational time on a laptop equipped with an Intel i7 processor (two cores, each one working at a frequency of 3.50 GHz) and 16 GB RAM. This same example was included in [17], achieving an eight = times lower computational cost. This means that the current version of the model allows simulating one hour of a real fire, like the one in the example, in one minute on a standard computer, which is highly competitive. In fact, although the computational cost depends on the size of the simulation area, the option of using different resolution levels enables the model to maintain a highly competitive computational cost in all types of simulations.

The three model parameter values for this simulation were as follows: mean absorption coefficient $a=0.095$ m${}^{-1}$, natural convection coefficient $H=15J{s}^{-1}$m${}^{-2}{K}^{-1}$, and correction factor of convective term $\beta =0.015$. The values of the fuel-type-dependent input variables used for the simulation are summarized in

Table 4. Fuel-type-dependent input variable values.

Fuel Type (NFFL/[34])	${\mathbf{M}}_0$	${\mathbf{M}}_{\mathbf{v}}$	${\mathbf{T}}_{\mathbf{f}}$	${\mathbf{T}}_{\mathbf{p}}$	${\mathbf{t}}_{1/2}$	C	${\mathbf{F}}_{\mathbf{H}}$	${\mathbf{F}}_{\mathbf{v}}$	${\mathbf{F}}_{\mathbf{s}}$
Timber grass (2/Pa-06)	$1.0$	$10\%$	1300	500	100	2000	$1.1100$	$0.4712$	$0.6759$
Brush (5/Eu-06)	$2.3$	$10\%$	1300	500	200	2300	$3.7780$	$0.5075$	$2.8280$
Dormant brush (6/Cl-02)	$2.2$	$10\%$	1300	500	200	2300	$3.3240$	$0.4888$	$2.6880$
Inflammable brush (7/Ea-08)	$2.4$	$15\%$	1300	500	300	2300	$3.9320$	$0.6752$	$3.0150$

for each one of the five fuel types in the simulation area, adapting the information from [34] to the NFFL classification [61]. The values of the maximum initial fuel load $M_0$ were obtained from [34]; the fuel moisture content $M_v$ depends on the fuel type and ambient humidity reported in [67]. The maximum flame temperature $T_{f,max}$ and pyrolysis temperature $T_p$ are standard values. The combustion half-life $t_{1/2}$ and heat capacity C varied slightly depending on fuel type, but had little influence on the ROS. Finally, the flame length factors, $F_H$, $F_v$, and $F_s$, were adjusted in the sense of least-squares to the data from [34] for the fuel types selected to represent the standard fuels according to the NFFL classification.

Table 5. Fuel-type-dependent input variable values.

Local Time	Temperature	Humidity	Wind Speed (m/s)	Wind Direction
3.45–5.00 p.m. average	32.02	27.17	2	260
5.00–6.00 p.m. average	32.01	27.43	3.3	300
6.00–7.00 p.m. average	31.50	27.43	4.75	300
7.00–8.00 p.m. average	31.42	27.50	4.75	360

summarizes the average hourly ambient temperature (Celsius), relative air humidity (%), wind speed (m/s) and direction (degrees from north), and wind gusts (m/s).

We compared the simulation of the fire spread considering the fire breaks shown in Figure 2 or not, and the results showed that these firewalls were key for saving the municipality of Osoño from the flames. Figure 3 and Figure 4 represent the burned area each hour and the fire front at 20.00, in both cases.

6. Discussion and Conclusions

We understand that research in applied mathematics is intricately linked to real-world issues. Interesting mathematical problems arise in engineering and other basic sciences, which pose challenges for mathematical research. Conversely, mathematics becomes an essential tool for improving our knowledge and understanding of complex systems. In particular, environmental issues constitute some of the major challenges for mathematical modeling because of their complexity and high degree of uncertainty. This is the underlying purpose of this work, without forgetting that technological advances allow faster progress in the development and resolution of complex models.

This article described the evolution of a mathematical model and the numerical resolution of a highly complex environmental problem involving wildfire spread. We studied PhyFire, from its initial concepts through to its current GIS-integrated version, paying special attention to the underlying mathematical, numerical, and computational aspects of this ongoing research work.

Throughout the entire process, we sought to strike a balance between model accuracy and computational cost, as the ultimate goal has always been to provide a useful tool for wildfire management.

Regarding the mathematical model, we have been making step-by-step improvements in pursuit of a simple model that is capable of representing the most-relevant phenomena to adequately represent the main features of a wildfire. Two of the key simplifications involved, first, considering a 2D model, while maintaining certain 3D phenomena, and second, considering a one-phase model, namely the solid phase, by parameterizing the gas phase. We paid particular attention to the simulation of heat transmission by convection and radiation, designing various proposals for the latter. We have always been mindful of the factors that most influence the spread of a fire, such as weather conditions, especially wind, the slope of the land, and fuel (type, distribution, and moisture content), all without increasing the model’s complexity while applying original mathematical concepts such as the use of a multivocal operator to simulate the effect of fuel moisture content. We added sundry sub-models to improve the simulation of certain key aspects, such as flame temperature and flame height. We also added the random phenomenon of fire-spotting, which is extremely dangerous during a real wildfire.

Concerning the numerical solution of the different models proposed, the objective was to maintain the lowest possible operational cost without compromising efficiency. We used various numerical methods ranging from finite differences to classical, mixed, and adaptive finite-elements. We also adapted the typical schemes of other types of problems to our models, such as duality methods. After checking different numerical schemes, the one that provided the best results in the current model was based on the use of P1 finite-element approximation on a regular mesh for spatial discretization and a predictor–corrector finite difference scheme for time discretization. In addition, certain ideas such as Actives Nodes, Radiation Molecule, and the ability to use several resolution levels were incorporated to maintain a fire’s simulation time well below its real-time evolution. This makes our simplified physical model highly competitive compared to other experimental models, and it can be adapted to simulate examples of different magnitudes.

From a computational perspective, the final code was optimized, and parallel computing was considered. Significant improvements were also made to the pre- and post-processing phases by integrating GDAL/OGR utilities into the code. These two steps are essential for completing the model and making it a fully functional and useful tool for wildfire management services. This was achieved by integrating the model into GIS and providing the necessary geospatial data.

Finally, the simulation results of the real fire reported here were highly accurate with a very competitive computational cost, whereby PhyFire was proven to be an effective tool for wildfire simulation.

This work is continuously evolving, and we are constantly working on improvements and solving new challenges. We are currently analyzing the effect of FMC on the simulation results through the multivalued operator. The positive results suggest that the use of remote sensing for providing reliable FMC data significantly improves the model´s efficiency. We intend to address other improvements in the physical model, such as the effect of soot on radiation. We also plan to address the complex issue of adjusting the model’s parameters, paying particular attention to the latest fuel classification systems and mapping available in Europe. This major task involves the simulation of a wide and varied set of real wildfires, for which all the necessary information is being gathered. We are also working on improving the coupling of PhyFire and HDWind, as we considered that studying the effect of fire temperature on local winds may be highly expedient. Lastly, a complex challenge involves coupling PhyFire with an atmospheric dispersion model to simulate wildfire smoke dispersion.