Effects of wind and tree density on forest fire patterns in a mixed-tree species forest

ABSTRACT It is well known that global climate change causes an increase in forest fire frequency and severity. Thus, understanding fire dynamics is necessary to comprehend the mitigation of the negative effects of forest fires. Our objective was to inform how fire spreads in a simulated two-species forest with varying wind strengths. The forest in this study was comprised of two different tree species with varying probabilities of transferring fire that was randomly distributed in space at densities (Ctot) ranging from 0.0 (low) to 1.0 (high). We studied the distribution pattern of burnt trees by using local rules of the two-dimensional model. This model incorporated wind blowing from south to north with strength (Pw) ranging from 0.0 (low) to 1.0 (high). Simulation results showed that when Ctot > 0.45 the fire covered the entire forest, but when Ctot ≤ 0.45 the fire did not spread. The wind effect on the variation of the amount of the burnt tree was maximized at the critical density and dramatically decreased with increasing Ctot. Additionally, we found that the term of Ctot and Pw plays an important role in determining the distribution.


Introduction
Biomass burning, in which both living and dead vegetation are burned, can be a process that emits greenhouse gases, reactive gases, and aerosols (e.g. CO 2 , CH 4 , and CO) into the atmosphere. This burning affects biogeochemical cycling of carbon, nitrogen, hydrogen, and oxygen between the biosphere and atmosphere, with global implications (Prather et al. 2012). In addition, individual forest fires that collectively contribute to biomass burning in some cases act as major disturbances to ecosystems and can negatively affect local economies (Mindas 1998). Given these negative impacts of fire, and the likelihood that the frequency and severity of fires are predicted to increase with global climate change (Matthew et al. 2010;Wotton et al. 2010), a better understanding of complex fire dynamics at the regional level using a multidisciplinary approach is imperative. In this study, our objective was to use simulation modeling to inform one component of fire dynamics -fire spreadin a simulated forest.
Spatial trends in the way fire spreads can result from complex interactions between weather, ignition, vegetation type, fuel moisture, and topography (Hargrove et al. 2000), making such patterns difficult to measure empirically. As a result, simulation studies are typically used to model how fires spread and to explain the patterns observed in the field. Such simulations are generally classified as either stochastic or deterministic (Halada and Weisenpacher 2005). Stochastic models are typically based on observational data and can incorporate variations in relevant burning parameters (e.g. fuel type, fuel moisture, wind) using statistical methodologies. These models have been used to explore the consequences of theoretical or assumed conditions on various processes. They have been validated against field or laboratory data, but often are not. Stochastic models require copious data on the interactions between burning parameters obtained from experimental fires, which are difficult to measure and strongly dependent on local environmental conditions. In contrast, deterministic models use global energy balance equations and are based on the assumption that energy transferred to unburned fuel is proportional to energy released by fuel combustion. These models take into account one or several processes of energy transfer from the burning zone to unburned fuel, which is helpful in understanding fire spread dynamics. The processes enable us to implement for planning purposes and for use in computer simulations to aid in prediction of the future behavior of existing and potential fires. In addition, the processes can provide the ways to incorporate variables inherent in the systems (e.g. animal trajectories, weather). In other words, this approach allows for the development of effective operational tools for forest fire simulations under real conditions. However, these models often lead to differential equation systems that require sophisticated and time-consuming numerical calculations and advanced high-performance computing environments (Boychuk et al. 2009). In addition, energy balance equations can be complicated in cases where many environmental factors (e.g. meteorological conditions, topography) must be incorporated to more accurately predict fire spread. Overall, it may be said that the difference between the two types of models is that the stochastic models can provide a range of potential outcomes with associated probabilities, whereas the deterministic models provide a single forecast based on a set of initial conditions and deterministic biophysical processes.
One alternative to deterministic and stochastic approaches discussed above is the use of cellular automata (CA) models. These models are computationally simple but capable of simulating complex behaviors (Berjak and Hearne 2002) and have been widely applied in ecological modeling (Matsinos and Troumbis 2002;Loibl and Toetzer 2003). Encinas et al. (2007) allowed for the propagation of fire between diagonal neighboring cells in CA models. Using the model, the authors predicted the spread of a fire front in homogeneous and inhomogeneous environments, and compared the results with those obtained from other simulation models. Yassemi et al. (2008) integrated geographic information systems (GIS) and CA techniques to develop a more accurate fire behavior model. The results from the simulation of the GIS-CA model and wave propagation spread models indicated comparable agreement. In addition, the authors validated their model of fire spread against empirical spatial data from the 2001 Dogrib Fire near Nordegg in Alberta, Canada to show that GIS-CA models can be used to accurately simulate real fire spread dynamics. Alexandridis et al. (2008) suggested a CA model describing the dynamics of a forest fire spread on a mountainous landscape. The model successfully explained the case of the wildfire that swept through Spetses Island in 1990. In addition to these studies, other models have been employed to model fire spread taking into account environmental variables (e.g. wind patterns) and heterogeneous landscapes (Li and Magill 2001).
Although CA models have been successfully employed to understand fire spread (Trunfio 2004;Russo et al. 2014;Ghisu et al. 2015), most models in the literature only consider a forest comprised of one tree species and either do not incorporate wind effects or do so through simplified wind metrics (Encinas et al. 2007). To our knowledge, few studies have developed CA models that consider both a multi-species forest and realistic wind dynamics.
In this study we simulated fire spread in a forest comprised of two tree species using a CA model that incorporated wind direction and strength. Overall, our model has an advantage when compared with those previously mentioned. In most existing models (both stochastic and deterministic), fire spreading has been simulated based on the (normal or abnormal) diffusion equation. Our model, however, mimics fire dynamics by considering the state of neighboring sites of a burning tree, giving advantages not only in computational cost but also in algorithmic efficiency. For example, modeling a simulated forest consisting of various physical properties of fuels, such as moisture, wind, and the geological structure, would require a set of complicated diffusion equations. However, the approach used in our model provides a simple way to contain the various properties regardless of the number of the properties. We used our model to understand how the distribution of burned trees varied depending on wind parameters and on the ratio and density of the two tree species in our simulated forest.

Model description
Our spatially explicit CA model depicted space in a square, two-dimensional lattice of size L £ L with periodic boundary conditions to minimize the boundary effect (Stevens 2000;Lee et al. 2006). The lattice is composed of cells where cell i, j represents a discrete position in the lattice at column i D 1, 2, …, L and row j D 1, 2,…, L. We used a Moore neighborhood that assumes that a neighborhood N(i, j) is defined by a central cell (i, j) and the eight cells (i Ã , j Ã ) surrounding it in the four cardinal and four diagonal directions such that (Kari 2005;Quartieri et al. 2010): (1) Each cell has a user-defined probability of changing its state according to its current state and that of its eight neighbors.
In this model we populated the lattice with a simulated forest comprised of two ʽspeciesʼ of trees that varied in their ability (or probability) of transferring fire (P t ): a ʽsusceptibleʼ tree species that had a high probability of transferring fire to and from other trees (P t D 1.0); and a ʽresistantʼ tree species with a low probability of transferring fire to and from other trees (P t D 0.5). In fact, real trees in a forest are likely to have various P t values. The values could also be changed according to the environmental conditions, such as moisture and temperature. The P t values in this study were selected to simplify the model algorithm and to make an easy analysis of the forest fire. Each cell within the lattice could be empty or populated by susceptible, resistant, burning, or burnt trees. In simulations, we varied the density of trees in the lattice (C tot ) in increments of 0.05 ranging from 0.05 (low density) to 1.0 (high density). The occupancy of a given cell n(i, j) was binary and determined according to: nði; jÞ D 1; when r ði; jÞ < C tot ðfor total treesÞ 0; otherwise (2) where r(i, j) was a randomly generated number at cell (i, j) by using the uniform probability function in [0, 1]. Here, a cell was occupied if its random number was less than the allowed total density of trees in the lattice for that simulation run. After the distribution, some of the trees were randomly chosen with a ratio C s (D [0, 1]) with respect to the whole trees. Thus, C s represents the density of the susceptible trees while (1-C s ) represents the density of the resistant trees. At the start of each simulation (time step t D 0), we introduced a single burning tree at the center of the lattice landscape as an ignition source (Berjak and Hearne 2002). At time step t D 1, fire could then spread from that cell to its neighbor cells containing a tree under no wind ( Figure 1A, B). The fire spreading was characterized as a probability (P t ). The value of P t indirectly reflects the total property of trees, such as tree type, density, and so on. In addition, any tree that was burning in one time step became a burnt tree in the following time step.
In our simulation, we assumed that wind impacts on the fire spreading. When the wind blew from south to north, trees in the three neighboring cells located above the center cell had a higher probability of transferring fire (P) depending on wind strength (P w ) ( Figure 1C). P w can be considered as the probability characterizing the degree of the impact. Consequently, the transferring probability,(P) can be simply written as: We varied wind strength from 0.0 (low strength) to 1.0 (high strength). We then simulated the model with all combinations of C tot , C s , and P w until no burning trees remained on the landscape, and we statistically averaged results over 100 independent model iterations.

Data analysis
To quantitatively measure the distribution patterns of burnt trees, we used three metrics: number of burnt trees divided by the system size, 200£200 (B); Moranʼs autocorrelation coefficient, often called Moranʼs I (I); and fractal dimension (F). I is an extension of the Pearson product-moment correlation coefficient to a univariate series and is a measure of spatial autocorrelation (Moran 1950;Fu et al. 2014) according to the following equation: where x i and x j are the x and y coordinates of the variable of interest, respectively. W ij is the weight between observation I and j, and S 0 is the sum of all values of w ij calculated as: Values for I can range from -1 (perfect dispersion) to C1 (perfect autocorrelation), and negative values indicate negative spatial autocorrelation while positive index values indicate positive autocorrelation. A zero value indicates a random spatial pattern.
Finally, F gives an index of complexity that compares how details in a pattern change with the scale at which they are measured (Pentland 1984). Although there are several different ways of measuring F, we employed the frequently used box-counting method (Wang et al. 2014) in which a square mesh of various sizes r is laid over the image containing the object. The number of mesh-boxes (N r ) containing a part of the image is counted, and fractal dimension F is given by the Figure 1. Graphical representation of the manner in which fire can spread from a burning tree in a central landscape cell (A) to trees in neighboring cells (B) according to our cellular automata model. The probability that fire will spread to a given tree (P t ) is based on the tree 'species' where one species is resistant to the transfer of fire (P t D 0.5) and the other is susceptible to the transfer of fire (P t D 1.0). The probability of fire transfer increases for both tree species when winds blow from south to north with strength P w (C).
To further examine the effects of wind on forest fire patterns, we calculated B, I, and F from results of simulations with varying combinations of C s and C tot and a fixed moderate value for wind strength. We plotted index values as functions of total tree density and measured two factors that characterized wind effects on the forest fire: (i) the increasing or decreasing slope in the fitting domain; and (ii) variation in B, I, and F for the different values of C s . For plots B vs. C tot , I vs. C tot , and F vs. C tot , we used fitting domains 0.45-1.0  (increasing region), 0.45-0.6 (decreasing region), and 0.0-0.6 (increasing region), respectively. The fitting domain for C tot was selected in consideration that when C tot < 0.45 (the critical density), fire could not spread at all. In the same reason, the fitting domain of I was set as [0.45, 0.6]. The end value of the domain (D 0.6) was selected because when C tot > 0.6 most trees were burned and the surviving trees were almost randomly scattered. In other words, for C tot > 0.6, Moranʼs I index does not contain the information regarding the spatial distribution of the burnt trees. Likewise, for C tot > 0.6, the values of F for the distribution of the burnt trees also has »2.0 because the burnt trees almost covered whole space, which means that no spatial structure of the burnt tree appeared in the domain.
We employed separate fitting functions for the slope of each index: log(aC tot ) for B, exp(-C tot /b) for I, and 1/(1Cexp (-C tot /g)) for F. Variation in index values was calculated by measuring the width of the distribution of B, I, and F at fixed values of C tot and was notated as H B , H I , and H F for variations in B, I, and F, respectively.

Results
Distribution patterns of the burnt trees for simulations with all combinations of C s (0.0, 0.5, 1.0) and P w (0.0, 0.5, 1.0) where C tot D 0.6 are shown in Figure 2. When wind strength and the density of susceptible trees were low (P w D 0.0; C s D 0.0), only small quantities of our theoretical forest burned, and they did so in distinct clusters. As the density of susceptible trees increased to a threshold of 0.5, individual trees were burned throughout the landscape, and the remaining unburned trees were clustered in lattice space. At the highest density of susceptible trees (C s D 1.0), most trees in the landscape were burned. The wind had the effect of increasing the size of burned areas as strength was increased from 0.0 to 0.5 to 1.0, and that effect was magnified with increasing density of susceptible trees. However, at the highest wind strength (P w D 1.0), distributional patterns of burnt areas were statistically the same regardless of the density of susceptible trees. Figure 3 shows metric values for B, I, and F for all combinations of P w (D 0.5), C s (0.0, 0.05, 0.1,…, 1.0), and C tot (0.0, 0.05, 0.1, …, 1.0). We found that B increased with higher values of C tot and P w , especially at the point where C tot D 0.45. This point represents a critical value where the landscape transitions from an unburnt forest to a burnt forest. Values for I also varied with levels of C tot (for C tot < 0.45). At higher levels of C tot (0.45 (C tot 0.6), I was relatively high regardless of wind strength, indicating that burnt trees were spatially clumped. At the highest levels of C tot (> 0.7), I decreased because burnt trees were located throughout the landscape and were, therefore, no longer clustered in distinct groupings. Finally, values for F slowly then dramatically increased with increasing C tot up to C tot > 0.6, where the entire landscape burned and F became saturated. Figure 4 shows the plots of B, I, and F against total tree density at varying densities of susceptible trees (C s D 0.0, 0.05, …, 1.0) assuming a fixed value for wind strength (P w D 0.5). We found that two slopes occurred in the plots of a against C tot : 1.580 §0.04 (r 2 D 0.9978) for P w 0.6 and 0.339 § 0.030 (r 2 D 0.9605) for P w > 0.6. Values of 1/b and 1/g increased linearly with increasing wind strength. Values for b and g were 13.412 (r 2 D 0.9912) and 9.205 (r 2 D 0.9976), respectively ( Figure 5). Finally, plots of H B , H I , and H F against total tree density assuming varying wind strengths are shown in Figure 6. When C tot 0.4, values of H B and H I were relatively constant because the fire could not spread beyond the small clusters of burnt trees. Values of H B , H I , and H F dramatically increased with the increasing of C tot with maxima at C tot D 0.45; however, at this threshold, maximum values of H B , H I , and H F decreased with increasing wind strength.

Conclusion
We used a CA model to understand the fire spread in a theoretical forest comprised of two tree species that varied in their capacities to transfer fire. We specifically used spatial statistics (i.e. B, I, and F) to evaluate changes in fire spread patterns with varying tree densities and wind strengths.
We determined that a total tree density of 0.45 represented a critical value; when tree density was above this level, fire was able to spread throughout the entire landscape unhindered. According to percolation theory (Janssen and Tauber 2005), a phase transition is predicted to occur when the tree density reaches this point, with the landscape moving from a relatively unburnt forest at tree densities below 0.45 (with low index values for F and B) to a relatively burnt forest at high tree densities (with high index values for F and B). In addition, we also found that variation in index values was maximized at a total tree density of 0.45, and the effects of wind were strongly dependent on the mixing ratio of the two tree species in our theoretical forest at this point. At tree densities above this value, the effect of wind strength on fire spread became low as closely spaced trees could easily spread fire even without strong winds to help the transfer. In particular, I values were maximized near 0.45, which means that the distribution of burnt trees forms a heterogeneous structure (including the spatial clusters of the burnt trees), not random. According to the definition of I, I value becomes zero for the case of no burnt trees (t D 0). When most trees were burnt (i.e. with only a few trees surviving), the value of I also goes to zero. Thus, it can be said that no information regarding the fire spreading is contained for too low and too high values of C tot . We found that distribution patterns of fire spread (represented by burnt trees), according to index values B, I, and F, could be explained with functions log (aC tot ), exp(-C tot /b), and 1/[1Cexp(-C tot /g)], respectively. These functions showed that fire spread patterns were largely governed by the product of total tree density and wind strength (C tot and P w ) and that the degree of fire spread was maximized at C tot D 0.5 and P w D 0.5. Thus, when total tree density is too low or too high, the wind effect has no impact Figure 5. Plots of a, b, and g (which represent parameters used in fitting functions) as a function of wind strength (P w ). Functions log(aC tot ), exp(-C tot /b), and 1/ (1Cexp(-C tot /g)) were used for fitting of spatial indices B (the number of the burnt trees divided by system size), I (Moran's I), and F (fractal dimension), respectively. on fire spread patterns. At low tree densities, the wind would need to be very strong to be able to transfer fire between widely spaced trees; however, at high densities, fire spreads easily from tree to tree even at low wind strengths. These patterns were further supported by the plots of variation in index values.
We purposely developed a simple model of fire spread to reduce the number of model assumptions, allowing for an analytically tractable simulation. In the model, we assumed that the wind speed and wind direction are constant over the entire simulation time period. However, in reality, wind blows in various directions with time-varying speed in a forest. In addition, our theoretical forest consisted of only two types of trees in order to simplify the model algorithms. When considering further ecological applications, this model could include additional parameters such as meteorology, vegetation, and topographical factors to improve model realism. For example, the effects of slope on fire spread could be represented by incorporating directional bias that approximates an anisotropic diffusive process across the lattice. Other species could also be added to the model with varying capabilities of transferring fire to further incorporate vegetation composition in models of fire spread. Although somewhat simplified, our model could have practical use in forest management to minimize fire risk.