The Forest Fire Model Revisited

TL;DR

This paper revisits the forest-fire model of self-organized criticality, showing through simulations that it exhibits behavior similar to rain and brain activity, suggesting relevance to real systems that hover near criticality.

Contribution

The study demonstrates that the Drossel-Schwable forest-fire model captures key features of systems near criticality, challenging the notion that exact scaling is necessary for relevance.

Findings

Model exhibits behavior similar to rain and brain activity.

Finite size scaling analysis performed on residence times.

Model's lack of exact scaling may still be physically relevant.

Abstract

Since Self-Organised Criticality (SOC) was introduced in the 1987 both the nature of the self-organisation and of the criticality remains controversial. Recent observations on rain precipitation and on brain activity suggest that real systems may dynamically wander about in the vicinity of criticality rather than tune to a critical point. We use computer simulations to study the Drossel-Schwable forest-fire model of SOC and find that it exhibits behaviour similar to that found for rain and brain activity. In particular we analyse the residence time for different densities of trees and perform finite size scaling analysis of the with of the distribution of residence times. We conclude that despite of the long known fact that this model does not exhibit exact scaling and power laws its behaviour may exactly for that be very relevant to real physical systems.

Tables

Table 1. Table 1: The total probability that the system is at a tree density more than one standard deviation from the average

L	Support
20000	0.327
15000	0.309
10000	0.320
5000	0.318

Full text

Self-Organized Criticality, power laws, temporal fluctuations, residence time, finite size scaling

I Introduction

Much of the research inspired by Self-organised Criticality took up the mantle from the paper by Bak, Tang and WiesenfeldBak et al. (1987) and studied how various simple dynamical systems may drive themselves into a critical state. Reviews of this report can be found in Jensen (1998); Pruessner (2011); Watkins et al. (2015). We are inspired to return to the discussion concerning the nature of the self tuning to a critical state or to the vicinity of such a state by the similarity found when analysing the size distribution of rain showers Peters and Neelin (2006) and the bursts of brain activity measured during fMRI scansTagliazucchi et al. (2012). Both studies find indications of critical behaviour in terms of approximate power laws and even features reminiscent of peaked, or perhaps diverging, fluctuations or susceptibilities. But both studies also investigated the time spent at different values of the control parameters (water vapour and number of voxels activated above threshold respectively). In both cases the distribution of residence times, the amount of time spent at a certain value of the control parameter, is found to exhibit a broad peak centred about what appears to be a critical value of the control. The immediate interpretation seems to be that even for systems as big as the atmosphere or with as many components as the brain the dynamics consists in some kind of feedback mechanism that is able to bring the system into the vicinity of a critical point. However, the dynamics couples the value of the control parameter to the fluctuations in a manner that moves the system around in this vicinity of the critical point. This is similar to suggestions previously put forward such as Zapperi et al. (1995); Sornette (1992). Intuitively, one may imagine something like the following in the case of precipitation: Nucleation of drops happens a critical value of the vapour content. The vapour in the atmosphere over the ocean gradually builds up towards that value and sometimes over shooting may even occur before nucleation is seeded. When precipitation events occur vapour is removed from the atmosphere and one observe oscillations between subcritical and supercritical region. For some reason the coupling between the driving (vapour formation) and the response (precipitation) produces so large fluctuations that a precise tuning to the critical value of the control parameter is excluded.

It was very early realised that one of the paradigmatic models of SOC, the Drossel-Schabl Forest Fire Model (FFM), doesn’t exhibit exact scaling as seen for ordinary equilibrium critical systemsGrassberger (1993); Pruessner and Jensen (2002); Grassberger (2002) and that this may be related to such feedback dynamics. Hitherto, most people have probably considered the lack of scaling in SOC models (the same complicated scaling is also seen for the original BTW sandpile model (See Chap. 4 in Pruessner (2011)) as a sign of broken promise in the sense that SOC suggested that the scenario of criticality and scaling seen in equilibrium systems were generic for a broad range of driven systems. Here we reanalyse the FFM with the observation on precipitation and brain activity in mind and we suggest that the model exhibits behaviour which may indeed be of relevance to the type of approximate, or near, critical dynamics seen in some real systems.

In this paper we simulate the FFM using the efficient algorithm described in Pruessner and Jensen (2002). We run the model in the usual mode where the the two probabilities with which empty sites become occupied by trees and the probability that a tree catches fire (from lightning) are prescribed and fixed. We then monitor the density of trees, the time the system spends at a given density, this we call the residence time, and the joint distribution of the sizes of tree clusters and the density of trees on the lattice. We analyse the behaviour as a function of the linear size $L$ of the lattice and find that the residence times to a very good approximation exhibit a Gaussian peak width a width of peak $\sigma(L)$ about a density $\rho^{*}(L)$. We find that for large values of $L$ the average tree density converges like $\rho(L)=\rho_{\infty}+bL^{-\alpha}$ towards the value $\rho_{\infty}=0.3976\pm 0.0001$, consistent with previously reported values Pruessner and Jensen (2002) with the $L$ dependence given by $\alpha=2.033\pm 0.037$ (These results are for $\theta=1000$, see below for details). We notice that tis density is far below the critical percolation density of about 0.6 for a square lattice. The width of the peak decays with the linear size of the system essentially as $\sigma(L)=16.82/L$ also in agreement with the analysis in Pruessner and Jensen (2002). The probability to find the system at a tree density more than one standard deviation away from the average is for all systems sizes about 32%, so very close to Gaussian behaviour.

It is of course difficult to know how to relate the lattice size of the FFM to the effective number of components of real systems such as the atmosphere or the brain, but if the FFM resembles such systems, the slow narrowing of the peak of residence times and the substantial fluctuations beyond one standard deviation for all system sizes suggest that SOC dynamics even for large systems generates an exploration of a significant range of control parameters rather than a self-tuning to a specific value or an infinitesimal region about such a value.

II Model details

We consider a two dimensional square lattice of linear extension $L$. In the language of the FFM a lattice site can be in one of three states: It may be empty, contain a tree or it may host a tree on fire. In the original version of the modelDrossel and Schwabl (1992) a random site would be chosen with uniform probability. If the site is empty it will become a tree site with probability $p$, if it contains a tree it will spontaneously be turned into a fire site with probability $f$ (that is stuck by a lightning) and if a tree site is nearest neighbour to a site on fire, the tree will become a fire with probability one. Fire sites are turned into empty sites with probability one. Here we implement this procedure by a slightly different, but very efficient, algorithm used in Pruessner and Jensen (2002); Clar et al. (1994); Grassberger (1993); Schenk et al. (2000). The main difference between the two algorithms is that now trees are grown in chunks of $\theta=\frac{p}{f}$ between two lightning attempts, this does not affecting the final statistics Schenk et al. (2000). Here we recall briefly the implementation of the algorithm in pseudo-code.

FOREVER $\{$

REPEAT $\theta$ TIMES $\{$

choose randomly a site $s$;

IF( $s$ is empty) THEN $\{$ $s$ becomes occupied $\}$

$\}$

choose randomly a site $s$;

IF( $s$ is occupied) THEN $\{$

collect statistics;

burn the whole cluster related to $s$ ;

$\}$

$\}$

The statistics is collected taking into account the overall density of trees, $\rho$, and the size of the burnt cluster $S$. An exception is made for the data represented in Fig. 1, where the biggest cluster has been removed before the computation of the mean and variance of all the clusters that are present in the system at each burning step.

III Results

We now discuss how the dynamics of the Forest Fire Model is related to the behaviour observed in rain and brain data. See Figs. 1 and 3 in Peters and Neelin (2006), Figs. 3 panel C, D, E and F in Tagliazucchi et al. (2012) and Fig. 10 in Scott et al. (2014).

We start by looking at the behaviour of the order parameter as function of the density of trees. In the study of rain by Peter and NeelinPeters and Neelin (2006) this corresponds to the precipitation rate as function of water vapour and in human brain studies by Chialvo and collaboratorsTagliazucchi et al. (2012) this corresponds to the normalised size of the the largest activity cluster as function of the total number of active voxels. The analysis of the human brain was later extended to mice, see Scott et al. (2014); Fagerholm et al. (2016), and a phenomenology very similar to the one reported for the human brain was found for mice brains. Inspired by percolation, two choices for the order parameter appear natural. One is the normalised size of the largest cluster of trees, i.e. exactly the quantity used by Chialvo and collaborators in Tagliazucchi et al. (2012) another, equivalent, is the average cluster size. We plot both these in Fig.1 in the form of heat-maps. In the same figure we also plot the variance of the cluster size distribution, this corresponds to the susceptibility in equilibrium statistical mechanics (see e.g. Ma (1985)) which diverges at a critical phase transition. The equivalent variance were studied in both the rain and the brain analysis.

To quantify to what extent the dynamics is able to drive the system into a critical state the rain and brain studies both analysed residence times. That is the distribution of times the system is found at a given value of the control parameter. A broad peak was found in the rain study, see Fig. 3 Peters and Neelin (2006) centred below the value of the water vapour at which the precipitation (order parameter) picks up. An even broader peak was found for the residence times in the brain study, see Fig. 3 panel E in Tagliazucchi et al. (2012) and Fig. 10 in Scott et al. (2014). In brain study the peak appears roughly to be centred in the vicinity where the order parameter starts to assume a significant non-zero value.

Fig. 2 contains the distribution of residence times and the distribution of the order parameter for $L=100$.

It is very demanding numerically to obtain good statistics for the joint distributions depicted in Fig. 1, hence the small systems size. To be able to do bigger systems we turn to one point measures below, but before that we want to point out that the distributions in Fig. 1 do suggest that the onset of something like critical behaviour occurs in a broad range.

Next we look at size scaling and are in particular interested in how the qualitative phenomenology of the FFM compares with the observations for rain and brain as we increase the system size and how the peak of the residence time distribution (left panel Fig. 2) changes as we increase the linear size of the system.

We find that the peak of the residence time is to a very good approximation Gaussian with a skewness essentially equal to zero.

The average tree density as function of system size approaches rapidly the its asymptotic value, as see in Fig. 3. The behaviour is well fitted by $\rho(L)=\rho_{\infty}+bL^{-a}$ and $\rho_{\infty}=0.3976\pm 0.0001$.

The standard deviation of the residence time distribution is shown as function of system size in Fig. 4 and we find that $\sigma(L)$ is well fitted by $\sigma(L)=16.82/L$. This behaviour is consistent with the suggestion in Pruessner and Jensen (2002) that the variance may decreases as $L^{2}$, corresponding to the number of trees on the lattice effectively being uncorrelated so that $\rho=\sum_{i=1}^{L^{2}}x_{i}$ can be see as a sum of uncorrelated variable $x_{i}$, where $x_{i}=1$ i the site $i$ is occupied by a tree and $x_{i}=0$ otherwise.

Table 1 lists the probability to find the system at a density one standard deviation away from the average for different system sizes. Gaussian behaviour with a probability for this deviation of about 32% is found for all system sizes. The Gaussian form is also directly evident from the plot of the distribution of residence times in Fig. 5.

In Fig. 6 we show the variance of the cluster size distribution for different system sizes $L=15,000$ and $20,000$. The peak of the distribution is in both cases located at tree densities above the value at which the residence time peaks. The peak in the residence time is located at $\langle\rho\rangle$ (see Fig. 3). In Fig. 7 we show the location of the peak in the variance $\rho_{var}^{peak}(L)$ as function of system size together with $\langle\rho\rangle(L)$. For all system sizes the peak in the variance is above the peak in the residence time. Unfortunately it is difficult to conclude how $\rho_{var}^{peak}(L)$ and $\langle\rho\rangle(L)$ behaves as $L$ goes to infinity.

The tree density at the location of the peak in the residence time, given by the average $\langle\rho\rangle$, is roughly located in the region where the order parameter picks up as indicated by Fig. 1. The pick up of the order parameter obviously occurs over a very broad range of values of tree densities and in principle this range will depend on the size of the system. To get high precision statistics for this is difficult but our simulations indicate that the order parameter picks up over an extended range of tree densities even for large system sizes. The broad pick up region and the same relative position of peaks is observed in the rain and brain studies, for rain see Fig. 1 and 3 in Peters and Neelin (2006) and for brain see Fig. 3 panel E in Tagliazucchi et al. (2012) and Fig. 10 in Scott et al. (2014).

In SOC studies observed tendencies to power law behaviour of event size distributions have been taken as a signature of ”critical” behaviour Jensen (1998); Pruessner (2011) ever since the original BTW paper Bak et al. (1987). The rain and brain study we are inspired by likewise exhibit power law like distribution of event sizes. This was reported for rain in Peters et al. (2002) and for the brain activity see Fig. 3 panel F in Tagliazucchi et al. (2012) and Fig. 9 in Scott et al. (2014). In the FFM it is, of course, the distribution of burned tree clusters which exhibits this power law like behaviour and has been studied extensively in the literature, see e.g. Pruessner (2011) and for careful studies of deviation from proper power law behaviour see Grassberger (1993); Pruessner and Jensen (2002); Grassberger (2002). We want to add a comment to this debate seen in the light of the studies of rain and brain activity. In both studies it is at least tacitly assumed and in fact made explicit in Tagliazucchi et al. (2012) that the control parameter fluctuates between low values at which the system can be seen as subcritical to high values where the system finds itself in a supercritical state. To understand how the FFM behaves in this respect we plot in Fig. 8 log-log plots of the distribution of sizes of tree clusters conditioned on specific values of the density of trees together with a plot of cluster sampled over all densities. The latter plot is the one considered previously in the literature and our graph for all tree densities exhibits the same kind of approximate power law behaviour as, e.g., seen in Pruessner and Jensen (2002).

We note that the cluster size distribution sampled for the tree density $\rho=0.394$ is lightly below the overall density (see Fig. 3) and the density $\rho=0.410$ is slightly above the overall density. These two conditioned cluster size distributions could be said to exhibit under and over critical behaviour in the sense that their tail is a little bit above and a little bit below the overall cluster size distribution. A behaviour reminiscent of the one reported for brain activity in Fig. 3 panel F in Tagliazucchi et al. (2012). And, remarkably, the cluster size distribution sampled for $\rho=0.402$ coincides with the overall distribution (see insert in Fig. 8). This behaviour might inspire one to think of the FFM driving itself around in a region of a critical point located at $\rho=0.402$ and while doing so mixing sub and super critical behaviour together.

To what extend this scenario describes exactly the dynamics and statistics of the FFM is not entirely clear. We do not know if, with increasing system size, scaling towards exact power law behaviour is happening for the cluster distribution conditioned for some appropriately value $\rho(L)$. To answer this question we would need to do extensive size scaling of the conditioned cluster size distributions, which, unfortunately, is numerically very demanding.

IV Discussion

It is very difficult numerically to obtain with certainty the residence time distribution, the order parameter and the variance of the cluster size distribution. One will expect that if the residence time distribution narrows down to a delta peak positioned at the density at which the variance diverges and the order parameter picks up, then the infinite FFM would exhibit a critical state like the one we know in equilibrium systems at a $2^{nd}$ order phase transition.

However, it may very well be that the dynamics for all finite systems sizes move the control parameter sufficiently around to complicate scaling. An extra complication consist in disentangling the scaling as the dissipation rate $f/p$ is taken to zero from the finite system size effectsPruessner and Jensen (2002).

Though we do suggest that the FFM might be a useful laboratory to investigate behaviour seen in real ”SOC” systems such as rain and brain, we also acknowledge that the origin of the broad cluster size distribution is still far from evident. The model self-organise to a density $\rho\simeq 0.4$ which is far below the critical density for percolation $\rho_{c}^{perc}\simeq 0.6$. So it is not at all obvious that the tendency to power laws, or scale free behaviour, see in the FFM is induced by mechanisms familiar from percolation. But given the similarity between the FFM phenomenology and the rain and brain analysis, understanding how approximative power law-like behaviour arises in the FFM, if not related to percolation, of course makes a better understanding more urgent.

**Acknowledgement

**LP gratefully acknowledges an EPSRC-Roth scholarship from the Department of Mathematics at Imperial College London. HJJ thanks Si Chen for interaction at the early stage of this project and is delighted to acknowledge a very helpful conversation with Dante Chialvo.

**Author contributions

**Both authors developed the conceptual ideas behind the research. L.P. performed the simulations. Both authors discussed the simulation results and jointly wrote the paper.

**Conflict of interests

**The authors declare no conflict of interest.