A Multi-Scale Network with Percolation Model to Describe the Spreading of Forest Fires

Abstract

Forest fires have been a major threat to forest ecosystems and its biodiversity, as well as the environment in general, particularly in the Mediterranean regions. To mitigate fire spreading, this study aims at finding a fire-break solution for territories prone to fire occurrence. To the effect, here follows a model to map and predict phase transitions in fire regimes (spanning fires vs. penetrating fires) based on terrain morphology. The structure consists of a 2-scale network using site percolation and SIR epidemiology rules in a cellular automata to model local fire Dynamics. The target area for the application is the region of Serra de Ossa in Portugal, due to its wildfire incidence. The study considers the cases for a Moore neighbourhood of warm cells of radius 1 and 2 and also considers a heterogeneous terrain with 3 classes of vegetation. Phase transitions are found for different combinations of fire risk for each of these classes and use these values to parametrize the resulting landscape network.

1. Introduction

Forest ecosystems and their biodiversity are currently threatened by forest fires, as well as the environment in general. In particular, in the Mediterranean regions, vegetation in rural landscapes is strongly subjected to the risk of wildfires, which can become a major cause of land degradation, with the consequent economic, human and patrimonial losses  [1,2,3,4].

To reduce the hazard of fire occurrence, understanding fire behaviour at several levels has been a core aspect of fire modelling throughout the years, with a special focus on the prediction of the fire spreading rate and the estimation of the released heat from the flame front [5]. Fire behaves according to three interacting physical factors: fuel availability (morphological and physiological characteristics of vegetation), weather (wind speed and direction, temperature, and relative humidity) and terrain (slope and aspect)  [6,7]—throughout this article these factors will be referred to as FWT conditions. These factors, along with the data of the initial fire condition, allow for the possibility to model fire behaviour [8]. Fire models such as Rothermel’s [9] predict a fire’s local behaviour, by using fuel model parameters as input. Fuel models, such as [10,11], are sets of parameters that describe the characteristics that classify certain groups of vegetation. Occasionally, these fuel models may include the environmental parameters of wind, slope, and expected moisture changes [9].

The field of percolation has played an important role in developing models and strategies to predict fire spreading. In parallel, cellular automata (CA) have been essential to overcome the complexity of natural dynamics of a given application area, as several studies show. For instance, work [12] presents a simple model which describes the propagation of a fire in a densely packed forest from a central burning tree, finding a critical probability value of approximately $1/3$ for the triangular lattice. Work [13] presents a bond percolation framework based on the medium ignitability and combustibility, showing the qualitative change of the rate of spread as a function of environmental factors and the fractal dimension of the burned area. The model uses CA with SIR states and von Newman neighbourhood, and includes heterogeneity as a random component. In work [14] a model is presented that predicts fire spreading in both homogeneous and heterogeneous forests, incorporating weather conditions and land topography. The algorithm used for the determination of fire fronts was found to be in agreement with real case scenarios.

Still, the application of these models on real case scenarios can be troublesome, given the natural complexity inherent to the number of factors involved. These factors, which are to be used as input parameters of a study, generally influence the pattern of land use/cover change. Besides, those considered for a specific study depend on the ease of access to data of the application region. Therefore, the selection of spatial and non-spatial factors depends on both the goals and the application area of the study [15].

One of the existing strategies for fire spreading mitigation is the implementation of fire-breaks, which are gaps in the available combustible material that prevents the fire from advancing. It is a critical structure in not only confining the burned area but also for delaying the spread of the fire [16,17], hence easing the combat. In [17], the authors use cellular automata modelling to identify efficient fuel break partitions for fire containment and study the efficiency of various centrality statistics. Here, a flat terrain and a single type of vegetation are taken into account. In [18], the same authors present a two-level approach where the dynamics of fire spread are modelled as a random Markov field process on a directed network. The cellular automata model used to parametrize the edge weights consider GIS meteorological and landscape information data.

Laboratory experiments can be planned to simulate a small-scale fire, but with limitations. An interesting study [19] compares theory results [20] with laboratory simulations, using randomly placed matchsticks with ignitable heads. The theory predicts that at critical percolation a fire front decelerates, whereas experiments indicate acceleration. Given this discrepancy, the author concludes that percolation theory models of forest-fire propagation that use simple site percolation are unlikely to be accurate. On the opposite side, large-scale experiments have their limitations in terms of reproducibility as well. Still, percolation modelling of fire-spread is of considerable importance since it describes the transition regime between a fire extinction (spanning fires) and an uncontrolled fire spread (penetrating fires) [12,19].

To overcome scale limitations, some examples using multi-scale and multilayer networks [21,22,23,24,25,26] and cellular automata with different approaches have come up with important findings. For instance, work [27] applies multiplex networks to model fire propagation, simulating a 3-layer of possible fire development: ground, surface, and crown, where each node of the multilayer represents a group of trees. At a larger (landscape) scale, work [28] presents a multi-scale structure, where the nodes are local land patches, each with their own spreading dynamics and, as such, present different times for the fire spread.

This work follows the line of research of the fields of percolation and complex networks. The main goal is to find a fire-break solution for the area of interest, at the landscape scale, based on fire behaviour according to the FWT conditions, at the local scale. The area of interest chosen for this study is the region of Serra de Ossa, in Portugal.

In this multi-scale approach, the constructed model consists of a two-scale network structure, where each node presents its own spreading dynamics, based on a quantitative local description of the land. For the particular case of Serra de Ossa, the landscape network has 10 nodes, resulting from the division of the area of interest in 10 land patches, according to different FWT conditions in each of them, measured at the local scale.

The local scale is defined as the range in which it is possible to define the limits of a land patch according to the distribution of a set of measurable features. The landscape scale is defined as the scale at which each patch of land is the element of study. At the landscape scale, local characteristics are not taken into account, but it matters to understand the connectivity of the nodes. For the sake of simplification, in this work the FWT conditions used to model fire behaviour at the local scale are simply the type of vegetation and its respective density within each patch.

Setting the area of each node as a cellular automaton, fire spreading is simulated using a percolation algorithm to find the threshold between spanning and penetrating fire regimes, based on different values of vegetation type and density.

For the case where only two classes of vegetation are accounted, i.e., occupied/empty cells of the cellular automaton, the results obtained for the percolation threshold in a square lattice with Moore neighbourhood, $p_c\approx 0.407$, are consistent with the literature [29,30,31]. The case for a neighbourhood of warm trees for each burning cell with different action radius is considered, changing the threshold value to $p_c\approx 0.725$, and $p_c\approx 0.375$, for first and second neighbours, respectively.

According to the work [32], fire hazard consists of the easy ignitability and propagation of the fire, together with its difficult extinction and it depends on morphological factors. Fire risk results from the combination of fire hazard with the burning probability. For the analysis of a heterogeneous scenario, the terrain is differentiated in 3 classes of vegetation, given the respective densities in each land patch. A value of relative influence is attributed to each of these classes. Since the fire risk has a more embracing definition, the influence level of a certain vegetation class is expressed by a measure of its respective fire risk. When the risk of each class is varied during spreading simulations, phase transitions that separate both spanning and penetrating fire regimes are found. The landscape network is then parametrized with combinations of the values of fire risk found for which such a phase transition occurs.

The main goal is to find an efficient fire-break solution that mitigates the spreading, to help civil forces fighting forest fires. Several models aim to fully thin the vegetation on certain areas, disregarding the consequences for the ecosystem. The strategy presented here only demands the containment of certain vegetation classes under a critical value of risk, aiming to maintain the regime of extinguishable fires.

Given the application geography of this model, the complexity of the problem is limited by restricting FWT conditions to those specific to that area. Still, the spatial extent to which this model can be applied strongly depends on previous land knowledge, which implies a demand for the inclusion of other areas of scientific and technological expertise.

As future work, after calibrating the spreading network dynamics, it is necessary to increase the resolution of the classification spectrum of autochthonous vegetation. Furthermore, when the statistical aspects are governed by the distribution of fuel, then a site percolation model would be most obvious; if the random aspects of propagation are the most important, then a bond percolation model is more appropriate [19], specially when the propagation is governed by vectorial influence from meteorological and orographic factors. The possibility for a mixed site-bond percolation model is considered, such that both fuel density and variation in fire propagation can be modelled, to enrich the existing model.

This work is organized as follows. After a list of the resources used in Section 2, a description of the structure of the model is presented in Section 3 and the scales of its application are defined. In Section 4, the procedure for computational simulations is presented, including the initial conditions of the problem. In Section 5, the graphical and numerical results obtained are shown, followed by the respective discussion, in Section 6, conclusions, in Section 7 and discuss further work in Section 8.

2. Materials and Methods

For the recognition and establishment of the external perimeter of the region of Serra de Ossa as the interest zone for this study, the ortophotos with 25 cm of resolution of Continental Portugal of 2018 were used as reference, obtained through https://cartografia.dgterritorio.gov.pt/wcs-inspire/ortos2018?service=wcs&request=getcapabilities, accessed on 30 November 2021 (WCS: Web Coverage Service), from https://snig.dgterritorio.gov.pt/, accessed on 30 November 2021 (DGT-SNIG: Direcção Geral do Território-Serviço Nacional de Informação Geográfica).

For the recognition of the morphology of the target zone, the information of the land use/cover of Portugal 2018 (COS18) was available at https://snig.dgterritorio.gov.pt/ (accessed on 30 November 2021) DGT-SNIG website, as well as the respective http://mapas.dgterritorio.pt/atom-dgt/pdf-cous/COS2018/ET-COS-2018_v1.pdf (accessed on 30 November 2021) technical specifications.

The software ArcMap from ArcGIS Desktop, version 10.8.1, was used to build the vector and raster files used in our study. The raster images were cut along the limits drawn and converted the result into 10 vector files (10 polygons) with information regarding the area, centroid and morphology of the vegetation in each polygon. The vector files were converted again into raster files and used as input in Python 3.8. Each respective polygon was converted into a matrix and treated as such throughout the algorithmic process and data analysis. The conversion consisted of a 1:1 mapping between the raster pixels and the matrix entries. CA simulations of the site percolation problem were run on the state matrix, which had non null entries according to the shape of the respective polygon. Vegetation class values were uniformly distributed throughout each of polygon’s area, according to the respective density values for each class. The initial ignition of each polygon took place at the centroid, with coordinates $(i_0,j_0)$.

Computation of network measures, as well as graphs and plots not related to cellular automaton simulations, were carried out in Mathematica 12.2, generously offered by SYMCOMP 2021, the 5th International Conference on Numerical and Symbolic Computation: Developments and Applications, at 25–26 March, Évora, Portugal, under the Young Research Awards.

3. Model

In this section the structure of the model is elaborated upon, developed at the different scales of application, as well as the respective fire dynamics in each of them.

The approach to the problem of fire spreading consists of a network structure applied to two scales of spatial range, the landscape scale and the local scale. A set of vertices and edges constitutes the network associated to each scale. Each vertex of the landscape network corresponds itself to a local network.

This structure is applied to a certain geographic area of interest for our study, Serra de Ossa. The area of interest is divided in land patches, according to its morphology, to form the vertices of the landscape network. The result of that land division. Each land patch is then divided into cells in a square lattice layout, to form the vertices of the local networks. Examples of internal cellular automata of some isolated patches can be found in Figure 1, Figure 2, Figure 3 and Figure 4. The whole set of cells of each land patch is ruled by the associated local dynamics.

The criteria used for patch division is based on the definition of local and landscape scales: The local scale corresponds to the spatial range until which it is possible to define a set of measurable terrain characteristics; the landscape scale is the scale at which each land patch, defined in the local scale, is the element of study and, as such, the mention characteristics no longer provide a relevant analysis. This means that it is possible to delimit the land patches in terms of a measurable characteristic of the terrain, such as vegetation type, dead fuel load, slope, etc., according to the available data. In this approach, only vegetation type an density were considered for patch definition.

A description of the terrain dynamics at each of the scales follows in the subsequent sections.

3.1. Local Dynamics

As defined previously, the local scale depends on the possibility to outline and measure some land features, whether it’s a river flow, soil humidity or shrub density. For land division, this study considers the type and density of vegetation and does not take meteorological and orographic factors into account. Also, the inputs are average measures of each vegetation type and density for data processing of local dynamics within each patch. Therefore, in the division process, the goal was to group land patches such that it is possible to identify the different types of vegetation and measure their respective density within the patch. As a natural outcome, irregular-shaped land patches emerge, with its own spreading dynamics, based on their respective FWT conditions which, in this case, only include morphology.

This strategy defines vertex layout of the landscape network, whereas the dynamics within each one dictates the strength of its connexions. Note that the delimited patches are disjoint and there is no space in the study area that does not belong to any patch. The “divide-to-rule” approach overcomes the impossibility to predict and monitor the consequences of local phenomena at a larger scale.

The subdivision of each patch into square cells is represented by a state matrix where CA simulations of spreading take place. The criteria for cell contamination is based on epidemiology model SIR, which make the correspondence of the states “susceptible” (S), “infected” (I), “recovered” (R) (analogously, “unburned”, “burning”, “burned”) to each cell $c_k$, namely, $s_k=1$, $s_k=-1$ and $s_k=0$, respectively.

Simulations consist of successive Markov processes. A Markov chain model is a stochastic process that analyses the probability that a state of a system at time $t_1$ changes to another state at the next time step $t_2$, regardless of the state of the system at previous time steps [15]. The memoryless property of this process is what provides stochasticity to the model and leads to emergence of interesting patterns. At each step of the process, the rules for changing the state $s_k$ of every cell $c_k$ along the process can be expressed as

$$\begin{array}{cc}\hfill & \left(a\right)p(s_l(t+{\tau }_1)=-1|s_k\left(t\right)=-1\wedge s_l\left(t\right)=1)=p_{kl},{\tau }_1\in \mathbb{N},p_{kl}\in [0,1],\hfill \\ \hfill & \left(b\right)c_l\in M_r\left(c_k\right),r\in \mathbb{N},\hfill \\ \hfill & \left(c\right)p(s_k(t+{\tau }_2)=0|s_k\left(t\right)=-1)=1,{\tau }_2\in \mathbb{N},\hfill \end{array}$$

where rule (a) defines the transition probabilities $p_{kl}$, which consists of the probability of a cell $c_l$ be contaminated at time step $t+{\tau }_1$ by a neighbouring cell $c_k$ which is burning at time t; rule (b) demands that cell $c_l$ belongs to a Moore neighbourhood of radius r of $c_k$ for any valid contamination. The Moore neighbourhood of $c_k$ is given by the set of coordinates $M_r\left(s_k\right)=\{(x,y+r),$ $(x,y-r),(x-r,y),(x-r,y),$ $(x-r,y-r),(x+r,y+r),$ $(x-r,y+r),(x+r,y-r)\}$; rule (c) indicates that a burning cell $c_k$ at time t is considered totally burnt (i.e., $s_k=-1$) at $t+{\tau }_2$ with probability 1.

Values ${\tau }_1$ and ${\tau }_2$ express different time measures. The former accounts for the time spent between the instant a cell $c_k$ is infected and the instant it has the potential to infect another; in other words it’s the latency time. After ${\tau }_1$ time steps, neighbouring cell $c_l$ is able to change its state from $s_l=1$ to $s_l=-1$. The value ${\tau }_1$ is intimately related to the heat generated by the burning cell $c_k$. Naturally, the thermal gradient is continuous, but for such a discrete model it is assumed that the heat front generated by $c_k$ immediately affects the whole cell $c_l$. Value ${\tau }_2$ can be associated to the convalescence time, i.e., the time spent between the states $s_k=-1$ and $s_k=0$, and it’s directly related to the quantity and flammability of material in $c_k$, regardless of its capacity to contaminate the neighbourhood. The way each patch is subdivided, creating it square cells, determines ${\tau }_1=1$. For simplification, it is assumed that ${\tau }_2=1$.

Assuming that each cell $c_k$ has the burning potential $P_k$ and the generated combustion heat $Q_k$ as intrinsic attributes, then each probability transition $p_{kl}$ must be a function of those attributes, which are related to terrain properties and thus, can be studied,

$$p_{kl}^p=f(Q_k^p,P_l^p),$$

(1)

where p is the patch where the process occurs. Facing a vectorial influence, whether meteorological or orographic, $p_{kl}$ from Equation (1) can be expressed as

$$p_{kl}^p=f(Q_k^p,P_l^p,{\overrightarrow{\psi }}_{kl}^p),$$

(2)

where ${\overrightarrow{\psi }}_{kl}^p$ corresponds to the vector of intensity $∥{\overrightarrow{\psi }}_{kl}^p∥$ and direction $\overrightarrow{kl}$, at patch p. Note that even though there is a cell indexation for $\overrightarrow{\psi }$, an average value can be considered for the whole patch, as convenient.

3.1.1. Classical Percolation

The classical case considers a homogeneous terrain, with one vegetation type, or one type of occupied cells. Each state matrix is previously prepared so that only the corresponding land patch area is active for simulations. The states $S:1,I:-1,R:0$ are distributed throughout the active area, such that $s_k=1$, and empty cells correspond to state $s_k=0$, as  shown in Figure 1.

The percolation process is executed according to the rules in Section 3.1. At $t=0$ the centroid cell is activated with state $s_k=-1$, corresponding to the initial ignition. For a starting cluster with only the initial ignition cell, the steps of the algorithm are (Algorithm 1):

Algorithm 1: Classical percolation.

For each $c_k$ with $s_k=-1$, evaluate $M_r\left(c_k\right)$, Contaminate $c_l\in M_r\left(c_k\right)$, changing $s_l=1$ to $s_l=-1$, Add $\#\left\{c_l\right|s_l=-1\}$ to the existing cluster, Repeat previous steps, with the updated cluster as input.

In the homogeneous case, the stochasticity is introduced merely with the initial state distribution, since there is only one type of vegetation. The pattern of a generated cluster at a density near the phase transition is shown in Figure 2.

3.1.2. Warm Neighbourhood

A case for a neighbourhood of warm cells was considered, representing pre-heated vegetation close enough to the fire front. A new condition is introduced to evaluate cell contamination at each step, which consists of generating a random number $u\in [0,1]$ and comparing it to each transition probability, redefined as

$$p_{kl}=\frac{\sum \left(c_l\right|s_l=1)}{\#M_r\left(c_k\right)},c_l\in M_r\left(c_k\right),$$

(3)

for neighbourhoods $M_r\left(c_k\right)$, with radius $r=1$ and $r=2$. In Equation (3), $\sum \left(c_l\right|s_l=1)$ is the sum of all susceptible cells in the neighbourhood of a burning one and $\#M_r\left(c_k\right)$ is the cardinality of the neighbourhood of $c_k$. This condition acts as a measure of burning potential $P_l$ given the available fuel near a burning cell. The more susceptible cells there are, the higher the burning probability. In this sense, sparse cells may accelerate the extinction of a fire, whereas cells with higher connectivity may lead the process towards an uncontrollable spread regime [12,19]. For an initial cluster composed of the initial ignition cell, the algorithm executes the following steps (Algorithm 2):

Algorithm 2: Percolation with warm neighbourhood.

For each $c_k$ with $s_k=-1$, evaluate $M_r\left(c_k\right)$, Evaluate $p_{kl}=\frac{\sum \left(c_l\right|s_l=1)}{\#M_r\left(c_k\right)}$, Contaminate $c_l\in M_r\left(c_k\right)$, if random $u\le p_{kl}$ and change $s_l=1$ to $s_l=-1$, Add $\#\left\{c_l\right|s_l=-1\}$ to the existing cluster, Repeat previous steps, with the updated cluster as input.

The pattern of the generated clusters at phase transition is similar to the homogeneous case, but critical density of occupied states varies, as shown in Figure 3.

3.1.3. Heterogeneous Terrain

For the heterogeneous case, 3 vegetation classes are considered, corresponding to 3 possible states, $\{0,1,2\}$. Each state matrix was previously set according to the densities of each class for each respective land patch, as shown in Figure 4.

This classification aims to introduce differentiation in land morphology, based on 2 assumptions. Let $h_i$ be the level of relative risk of the i-th class:

$$\begin{array}{cc}\hfill & \left(a\right)h_0\equiv 0,\hfill \\ \hfill & \left(b\right)h_1<h_2,\hfill \\ \hfill & \left(c\right)h_i\in [0,1],i=\{0,1,2\}.\hfill \end{array}$$

Assumption (a) states that there exists a vegetation class that is not susceptible to fire occurrence, meaning that the corresponding state $s_k=0$ of cells $c_k$ never changes throughout the spreading process; assumption (b) states that the relative risk of one of the two vegetation classes susceptible to fire is lower than the other, regardless of its value. Assumption (c) constrains the limits of the scale of relative risk. The conversion to an absolute scale is possible, given the knowledge of the absolute extreme value of a fire risk scale.

For the heterogeneous case with 3 vegetation classes, for an initial cluster of 1 initial central cell in combustion, the percolation process follows the steps (Algorithm 3):

Algorithm 3: Heterogeneous percolation.

For each $c_k$ with $s_k=-1$, evaluate $M_r\left(c_k\right)$, Evaluate $s_l$ of each $c_l\in M_r\left(c_k\right)$, Generate random $u\in [0,1]$ and compare $u\le h_i$, If previous step is true, then $s_l=-1$, Add $\#\left\{c_l\right|s_l=-1\}$ to the existing cluster, Repeat previous steps, with the updated cluster as input.

Similar patterns of clusters generated in the homogeneous and heterogeneous cases are observed, which is natural, given the similar transition rules. However, for the heterogeneous case the vegetation class densities are fixed according to information on the morphology of the study area. Instead, the relative risk of classes varies, $h_1$ and $h_2$, given $h_0\equiv 0$, and combinations of $(h_1,h_2)$ are found for critical percolation. This strategy is particular useful for choosing what vegetation class to constrain, in general, in a certain region.

3.2. Landscape Network

The development of a model for local dynamics in Section 3.1 aims to generate a set of values related to the internal patch dynamics that determine its connectivity among its adjacent ones. Thus, the landscape graph is parametrized based on a function of initial FWT conditions.

The parameters consist of the edge weights that connect each pair of neighbouring vertices $(p,q)$ of the graph G,

$$G=(V,E),V=\{1,\dots ,N\},E=(p,q)|(p,q)\in V^2,p\ne q,$$

(4)

which can be represented by the adjacency matrix,

$$M\left(G\right)=\left[\begin{array}{cccccccccc}0& 1& 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 1& 0& 1& 0& 0& 1& 1& 1& 0\\ 1& 1& 1& 0& 0& 0& 0& 1& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 1\\ 0& 0& 0& 0& 0& 0& 1& 1& 1& 1\\ 0& 0& 1& 0& 0& 1& 0& 1& 0& 0\\ 0& 0& 1& 1& 0& 1& 1& 0& 1& 1\\ 0& 0& 1& 1& 1& 1& 0& 1& 0& 1\\ 0& 0& 0& 0& 1& 1& 0& 1& 1& 0\end{array}\right],$$

(5)

which is not weighted and, therefore, is symmetric. The adjacency of $(p,q)$ is valid for any two patches with any common border segment. The edges are parametrized by values of risk of the i-th vegetation class,

$$H_{pq}=f\left(h_i\right),i=\{0,1,2\}.$$

(6)

In the absence of any vectorial influence, it is assumed that any neighbour $q=1,\dots ,k$ has equal risk of being infected, that is, $H_{p1}=...H_{pk}$, where k is the degree of node p. Thus, in the simplest case, based only on a morphology analysis,

$$\begin{array}{c}\hfill H_{pq}=H_p\Rightarrow H_p=f\left(h_i\right),i=\{0,1,2\},\end{array}$$

(7)

and the representative graph, G, of Serra de Ossa, is shown in Figure 5. In the particular case of this study, $H_p$ provides us the relative risk values $(h_1,h_2)$ for which a phase transition occurs in patch p, always considering a null risk for class 0, there is, $h_0=0$. Any neighbouring patch q is affected by fire percolation throughout p. In the presence of wind or slope, the neighbours q would have different probability of being infected, namely those that are aligned with the direction of the vectorial influence.

It is possible, however, to specify the form of $H_p$. For simplicity, the suggestion for this study, is a linear combination of values $h_i$, given the i-th class density values for each patch p, $d_{pi}$. For this particular application,

$$H_p=\sum _{i=0}^2{d}_{pi}{h}_i.$$

(8)

Equation (8) can be changed according to what is most convenient for the study, and detailed according to the goals of experts from more applied research fields. In addition, future work developments regarding the inclusion of meteorological and orographic factors should consider neighbour differentiation, that is, $H_{pq}$.

4. Application Area

As mentioned previously, the application of this model to a real case scenario is strongly dependent on the natural complexity of the region and ease of data aquisition [15]. The region of Serra de Ossa is prone to fire occurrence, reason why the volume of studies and collected data on the target area is significantly higher in relation to other areas. Besides, its substantial elevation suggests the possibility to add a well-behaved vectorial component to the current model, in future work. Thus, the used FWT conditions are a consequence of applying this model to the particular region of Serra de Ossa. In this section the preparation of the study area for the execution of computer simulations is explained.

4.1. Delimitation of the Study Area

Using visual recognition based on satellite imagery, the external limits of the region of Serra de Ossa were drawn, defining an hexagonal polygon, to be subdivided afterwards, according to the internal land morphology. The representative area is shown in Figure 6.

4.2. Internal Division and Morphology

Once the external perimeter is established, the morphology of the area was described within it, aiming to group the existing vegetation into relevant classes, with different levels of fire risk. For the effect, the information of vectorial files from COS18 was essential, along with the respective technical specifications and GIS tools to create a map with all categories of land occupation/use. Every class of COS18 within the territory was then regrouped into 3 relevant classes:

Class 0

May include urban and water zones, and also agricultural crops regularly watered; fire risk level $h_0=0$.

Class 1

May includes agro-forestry surfaces and bushes; fire risk level $h_1$.

Class 2

May include forests and dense shrubs; fire risk level $h_2$.

Based on this reclassification, the resulting division in 10 land patches is presented in Figure 7.

Table 2. Percentage of vegetation classes 0, 1 and 2 for each land patch.

Patch	Class 0	Class 1	Class 2
0	0.17	0.79	0.04
1	0.11	0.58	0.31
2	0.32	0.56	0.12
3	0.16	0.73	0.11
4	0.09	0.06	0.85
5	0.11	0.12	0.77
6	0.04	0.77	0.19
7	0.09	0.29	0.62
8	0.07	0.67	0.26
9	0.02	0.0	0.98

discriminates the percentage of each class for each division.

The values of Table 2 allow the distribution of the set of states $s_k=\{0,1,2,-1\}$ by every cell $c_k$, to construct heterogeneous state matrices. An example can be seen in Figure 4.

5. Results

In this section the results of simulations performed are shown, first under the classical framework (Figure 8), then considering a neighbourhood of warm cells (Figure 9 and Figure 10) and, lastly, for the heterogeneous scenario (Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15), as explained previously in Section 3.1. The respective phase transitions for each case are shown. Every point in each plot corresponds to the cluster size value, averaged over 100 generated clusters in total.

5.1. Classical Percolation

A clear phase transition is observed for the classical case (Figure 8), in a homogeneous terrain. At critical density $d_c\approx 0.41$, fire spreading ceases its extinction regime to the uncontrollable regime, where the generated cluster increases with the patch size.

The mean cluster size as a function of density of occupied sites follows the same pattern for each patch, given the universality of the phenomenon of percolation, and only differs after the phase transition due to the different patch areas.

5.2. Warm Neighbourhood

Still in a homogeneous scenario, considering a warm neighbourhood of cells implies, in fact, introducing a restriction to the infection rules: not only neighbouring cells must be available for burning, $s_k=1$, as the neighbourhood they belong to must also be sufficiently populated with other susceptible cells. In this sense, the community matters and influences the strength of the spreading dynamics. For an interaction with $r=1$, i.e., first neighbours, a phase transition is observed for critical density value of $dc\approx 0.725$ (Figure 9), which is naturally higher than the classical case, since the added restriction only hampers infections at each step. When the radius of interaction is extended to $r=2$, the influence of the number of neighbouring cells considered for the contamination at each step superimposes the community-based restriction. Thus, the result is a lower value for the phase transition, at $d_c\approx 0.38$ (Figure 10).

5.3. Heterogeneous Terrain

For the heterogeneous case, 3 classes of vegetation were considered for the initial state matrix to perform the simulations. The respective density values for each patch were fixed according to Table 2, and the fire risk values varied for each class. Class 0 has a risk $h_0=0$ of being burned, considering that is related to urbanized or watery surfaces. The details for the classification of classes 1 and 2 is beyond the scope of this work and should involve expertise in other scientific research fields. It is only assumed that there is a class of intermediate fire risk $h_1$ and a class with the highest fire risk, $h_2$, such that $h_1<h_2$, hence, the horizontal axis in plots (Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15) presents a $h_2$ variation from $h_1$ upwards, in a normalized scale, with increments of $\Delta d=0.01$. A 3D graphical visualization is more accurate, but the choice of a 2D analysis allows a higher focus of phase transition patterns for specific $h_1$ values. As a consequence, the presented plots reveal more than one possible risk combinations. For patches where that occurs, the combination with lowest $h_1$ or for which the phase transition has a more significant pattern were chosen to be listed in

Table 3. $H_{pq}$ parametrization for phase transitions of graph G. Values $(h_1,h_2)$ are approximated.

p	$({\mathit{h}}_1,{\mathit{h}}_2)$	${\mathit{H}}_{\mathit{p}}={\sum }_{\mathit{i}=0}^2{\mathit{d}}_{\mathbf{pi}}{\mathit{h}}_{\mathit{i}}$
0	$(0.3,1.0)$	$0.277$
1	$(0.2,0.55)$	$0.2865$
2	$(0.4,0.60)$	$0.296$
3	$(0.3,0.50)$	$0.274$
4	$(0.1,0.29)$	$0.2525$
5	$(0.1,0.31)$	$0.2507$
6	$(0.2,0.70)$	$0.287$
7	$(0.1,0.41)$	$0.2832$
8	$(0.2,0.57)$	$0.2822$
9	$(0.1,0.25)$	$0.245$

.

Phase transitions in Figure 11 occur for patches 4, 5, 7 and 9 at $h_2$ values of approximately $0.29$, $0.31$, $0.41$ and $0.25$, respectively. Patch 1 starts a phase transition at $h_2\approx 1.0$, but the spreading never reaches the uncontrollable regime. For every patch, cluster sizes approach the patch vegetation density at the post-critical regime.

For $h_1=0.2$, a phase transitions occurs for patches 1, 4, 5, 6, 7, 8 and 9, respectively for $h_2$ values of $0.55$, $0.3$, $0.32$, $0.7$, $0.33$, $0.57$ and $0.25$. The post-critical regime show in Figure 12 for patches 1, 6 and 8 show that the growth of the mean cluster size is stabilized by the density of available fuel.

Figure 13 presents the results for $h_1=0.3$ and shows a phase transition for patch 3 occurring at $h_2\approx 0.5$, after which the mean cluster size stabilizes at approximately half of the patch area. Patch 0 starts a very smooth phase transition at $h_2\approx 1.0$, which is listed in Table 3.

In Figure 14, phase transition for patch 2 is not as abrupt as for the situation in Figure 11 or Figure 12, and also fire spreads throughout half the patch. Finally, with a value of $h_1=0.5$, all patches have experienced phase transitions and are in the regime of penetrating fire, as shown in Figure 15.

Table 3 shows the values for tuples $(h_1,h_2)$ with lowest $h_1$ and the results of respective linear combinations, from Equation (8), using density values $d_{p}i$ from Table 2.

With the parametrization resulting for the linear combination chosen for $H_{pq}$, some centrality measures are calculated, as shown in

Table 4. Some centrality measures for graph G (approximate values).

Patch	Eigenvector Centrality	Closeness Centrality	Betweenness Centrality
0	$0.070$	$0.5$	$0.0$
1	$0.070$	$0.5$	$0.0$
2	$0.135$	$0.75$	$16.467$
3	$0.120$	$0.692$	$9.6$
4	$0.081$	$0.563$	$1.467$
5	$0.144$	$0.75$	$10.733$
6	$0.138$	$0.75$	$18.467$
7	$0.050$	$0.474$	$0.0$
8	$0.093$	$0.563$	$2.8$
9	$0.098$	$0.563$	$2.467$

.

6. Discussion

The classical percolation algorithm provides a result that is in agreement with the literature and provides a good reference to compare any changes to the base model, such as considering warm neighbourhoods and 3-class heterogeneous terrain morphology.

In the homogeneous cases, it’s observed that after the phase transitions, the mean cluster size increases with the size of the patch. In a real case scenario, excluding other factors, mitigating fire spreading risk in such conditions would consist of containing the percentage of area covered with vegetation below the critical density value $d_c$. That would explain the low propensity for fire occurrence in land with a large extension of consistently watery crops (assuming a zero burning risk for these surfaces).

Accounting for a warm neighbourhood, the critical density has a higher value than that observed for the classical model where only first neighbour interactions are allowed, the same value is significantly lower when the radius is extended to second degree neighbours. In both cases, the pre-heating depends on the amount of available fuel, which depends on the number of susceptible cells in the neighbourhood. In reference to the classical framework, the warm neighbourhood scenario only adds up the evaluation for available burning material, which hampers propagation at each step. However, when second degree neighbours are considered, the number of cells superimposes that restriction, thus lowering the value $d_c$, easing propagation at each step. In a real case scenario, excluding other factors, this would mean that long-range interactions, such as trees capable of ember projection, for instance, would be determinant in maintaining conditions for uncontrollable fire spread. Given the non-linearity of forest fire phenomena, the study of community-based dynamics is essential to understanding the spread.

In the heterogeneous case, fire risk values are associated for each of the 3 classes of vegetation, with land patch densities listed in Table 2. The values found to parametrize the landscape graph (Figure 5) are listed in Table 3. Instead of assuming a binary existence of vegetation for each cell, it is only assumed that, regardless of the amount of vegetation, the fire risk associated falls into one of the class categories. Watery, rocky or urbanized surfaces, for instance, belong to class 0, with $h_0=0$ and, probably, terrain with some sparse low grass. Likewise, more complex forms of vegetation, well-classified by experts, most likely will belong to classes 1 and 2. Anyway, for these simulations, initial state matrices were considered to be totally covered by these three classes, in proportions respective to each patch. For known distributions of these classes, the variation of fire risk for classes 1 and 2 not only opens the opportunity for a more specific classification, depending on the application geography, as it allows the visualization of a pruning strategy for each class, as a fire spreading mitigation method. For instance, if a class is determinant in setting a phase transition at a certain risk level, a controlled pruning of the corresponding vegetation may be worth considering.

The application geography of this model, Serra de Ossa, is of particular importance, given its special propensity for fire occurrence. Thus, providing a fire risk-based network which allows mapping of such a region is a preventative approach to the problem of fire spreading that might be useful as a complementary effort to fight forest fires.

7. Conclusions

Based on our results, a fire-break solution could be approached by creating a pruning plan for vegetation of a certain class, to lower its associated fire risk. That would overcome the indiscriminate elimination of large extensions of vegetations, threatening forest ecosystems. Adequate measures would imply the maintenance of the spanning fire regime, coupled with a model of fire occurrence probability, as a preventative combat strategy.

8. Further Work

This work is part of a continuous study that aims to develop a preventative strategy to fight wildfires in the region of Serra de Ossa.

Given the natural complexity of the real case dynamics, several minor studies must be developed within this major project, namely the effect of land division on phase transitions. The division performed in this work aims at grouping patches in a way that is possible to define the type of vegetation and measure the density for each type. However, maintaining that criterion, changing the patch layout can change the results of phase transitions, whether in regular or irregular lattices. One of the goals is to find a suitable layout to study the landscape.

Another important study that must be developed is the application of the current model accounting the wind and slope to differentiate the preferred fire path throughout the patch. In that case, values $H_{pq}$ will then depend on neighbour p. The meteorological and orographic factors depend also on the vegetation structure, and first it is necessary to guarantee access to that respective data.

Finally, to calibrate the model it is necessary to compare the pattern of burned areas from previous real fire occurrences to the pattern obtained from generated clusters. Under the same FWT conditions, measuring the cluster size is important, but also measuring the fractal dimension of both perimeters, especially in local segments of the perimeter, with different incidence angles of the vectorial influences, and different FWT conditions, in general.