Electronic Supplemental Material for ‘Self-organised criticality and the response of wildland fires to climate change’

SALVADOR PUEYO

Departament d’Ecologia, Universitat de Barcelona, Avgda. Diagonal 645, 08028 Barcelona, Catalonia, Spain.

Tel. (34)934021506

Fax (34)934111438

E-mail:

S1. Contents of the Electronic Supplemental Material

The following supplemental material is included:

  • Table S-I, with the whole set of symbols and acronyms used in the paper and in the Electronic Supplemental Material. This is the extended version of Table I, which only contemplates the symbols used in the paper.
  • Table S-II, with the parameters used for the simulations.
  • The extended version of the materials and methods for Section 4 (‘Fire dynamics in model and in reality’).
  • The analytical and numerical work in which Section 5 (‘Analysis of the model using phase transition theory’) is based, including both methods and results.
  • The materials and methods for Section 6 (‘Response of boreal forest fires to climate change: a forecasting exercise’).
  • Table S-III, which belongs to the results of Section 6 (‘Response of boreal forest fires to climate change: a forecasting exercise’).
  • The extended version of Section 7 (‘The contrasting case of tropical rainforests’).
  • The figures, tables and references accompanying the above texts.

S2. Extended version of the materials and methods for Section 4 (‘Fire dynamics in model and in reality’)

I analysed the dynamics of the SOCFUS model and empirical data from wildland fires in the boreal forest dominated Forest Protection Area of Alberta (Canada).

I carried out the two simulation experiments in Section 3.3, using the algorithm in Section 3.2 and the parameters in Table S-II.

Alberta’s historical fire register is available in the Fire Database of the Government of Alberta, Land and Forest Service, Forest Protection Division ( This register covers 1966 to 2003. However, pre-1983 data are unreliable and biased toward large fires, so I only used the dataset from 1983 to 2003. All treatments involved the sizes of the fires, expressed in hectares. An examination of the dataset clearly reveals that data referring to fires smaller than 1 ha are completely unreliable, largely due to rounding effects, so I established a lower cutoff at s0 = 1ha (see Section 2.4). This left 5,388 fires for the analyses. There are also some problems caused by rounding for larger fires, but these were not so serious. Rounding can be detected because there are well defined peaks of frequency in the size distribution, corresponding to round numbers such as 10, 15, 20, 100, 150, etc. However, this is not so apparent when the data are binned as are the figures in this paper.

I represented the empirical probability density functions (p.d.f.) as follows: I divided the data into bins i satisfying (in hectares), for the integer . The (logarithmically) central size in a bin is ha, and its estimated probability density is , where ni is the number of fires in the bin, N is the total number of fires, and 2i is the width of the bin. In the figures I display log10(si) vs. . It is convenient to use multiplicative intervals like these for power laws, as the distribution of data into different bins becomes much more balanced than in common histograms.

I did not use maximum likelihood (m.l.e.) to estimate the parameter  of the power law distribution, as this method is sensitive to the aforementioned rounding effects and to any other departure from a strict power law. Such departures are small but meaningful; I deal with them in detail in Section 6. Moreover, m.l.e. gives the same weight to all fires, which is a problem in our case because: (i) most fires are very small, so small fires have the most determining effect on the result, (ii) most of the burned area is due to large fires, so it is not convenient to give them such a small weight in the estimation process, and (iii) small fires are the most severely affected by rounding problems. When empirically fitting a scaling function, it is preferable to use a method that gives equal weight to all scales and that is insensitive to small departures from scaling. I estimated  by nonparametric regression (Kendall’s robust line-fit method, Sokal and Rohlf, 1995), as the slope of the log-log representation of the p.d.f. is –,.

The few largest bins in the p.d.f. are the most severely affected by sampling noise, due to lack of data, and are often empty. This appeared to be particularly problematic for the sequence of regressions in Figure 5d-f, so, in this case, I just took the series from the first to the last bin before the first empty one, which meant neglecting 0-2 nonempty bins in each regression, out of 9-14 bins. I encountered the opposite problem in the sequence resulting from the simulations in Figure 5a-c: in this case, the data are so abundant that even the bending function is well sampled and produces a clear deviation from a straight line at the upper end of the distribution. I also eliminated the upper end in this case, though in a less systematic manner. The power laws in Figure 3 and in Section 6 were not subject to any special treatment before undertaking the regressions.

I took the size of the largest fire in each set as an estimator of the parameter sM of the power law. This is the m.l.e., assuming that Equation (3) holds, and it is robust to the aforementioned problems, but biased. It is appropriate for the objective of this section, which is the qualitative comparison of model and reality, but I develop and apply a more sophisticated method to obtain quantitative predictions in Section S5.2.3.3 (Methods for Section 6).

The power law is one of the three main distributions used to fit the asymptotic range of probability functions. The other two are the exponential and the bounded distribution (Pickands III, 1975). Most text-book distributions are asymptotically exponential (Pickands III, 1975). Therefore, I tested whether the p.d.f. of the set of fires in Alberta is better fitted by a power law or by an exponential function, with criteria similar to those I used to fit the parameters. The representation gives a straight line for a power law and a curve for an exponential, while gives a straight line for an exponential and a curve for a power law. I tested the exponential null hypothesis against the power law alternative as follows: (i) I calculated the squared Pearson’s correlation R2 for ; (ii) I fitted a straight line to (the use of either Pearson’s or Kendall’s regression did not affect the result of the test); (iii) I randomised the residuals 999,999 times, (iv) for each pseudoreplica, I calculated R2 for ; and (v) I obtained the significance level  that would allow rejecting the exponential distribution, which is the proportion of cases (pseudoreplicas plus empirical set) in which R2 is equal or larger to that in the empirical set. I carried out the symmetric procedure for testing the power law null hypothesis against the exponential alternative.

In addition to the overall statistical distribution of fire sizes, I analysed the response to weather fluctuations, in both the model and reality. In the case of the model, I took the data from experiment 1 (Section 3.3), which corresponds to a landscape subject to short-term weather fluctuations. In this experiment, every fire event develops under a given value of the environmental variable re, which is chosen at random each time lightning strikes and, if this causes a fire, is maintained until the fire goes out. I separated the fires into groups corresponding to regular intervals of re and analysed the p.d.f. separately for each group.

In the case of the empirical data, I took the Canadian Fire Weather Index (FWI) as a rough equivalent of re. This index is calculated exclusively from meteorological variables and is assumed to capture the effects of weather on fire intensity, mainly through its presumed consequences for fuel moisture, and also taking into account the direct effect of wind on fire (Van Wagner, 1987; Flannigan and Van Wagner, 1991). Until 1995, the Forest Protection Division of Alberta registered FWI values calculated from measurements taken in the closest meteorological station at the moment in which each fire was thought to have started. I took the fire data from 1983 to 1995, distributed the data in groups according to FWI, and analysed the fire size p.d.f. for each group. I did not take regular intervals of FWI because this would have produced quite an unbalanced distribution, preventing the efficient use of our limited data set. Instead, I arranged the fires into groups with a similar number of data per group. The ranges of FWI in each group are {1-3, 4-6, 7-9, 10-11, 12-13, 14-15, 16-18, 19-21, 22-24, 25-29, 30-37, 38-86}. I labelled each group with its mean FWI.

In the case of the model, I compared the effect of short-term meteorological fluctuations to the effect of sustained climate conditions, i.e. compared the outcomes of experiments 1 and 2. I also distributed the data in intervals of re and calculated the mean fire size in each bin. Unfortunately, it was not possible to subject the empirical data to a similar treatment. The time series is relatively short, and subject to significant modifications that are not related to either climate or forest self-organization. As I mention above, the data from 1966 to 1982 are of poor quality and are biased to large fires. The improvement in data quality in 1983 occurred at the same time as fire suppression was strengthened, in response to the large fires of the previous years (Nash et al., 1999). However, the effectiveness of fire suppression seems to have gradually declined between 1993 and 1998, due to a policy of decreasing government expenditure (Nash et al., 1999). These mid-term changes largely cancel each other out when the average response to short-term meteorological fluctuations is examined, but they obscure any long-term response.

S3. Methods for Section 5 (‘Analysis of the model using phase transition theory’)

A second-order phase transition cannot strictly be observed in systems with a finite size, because they cannot display an infinite fire. However, we can test whether Equation (5) is satisfied while fires are still reasonably smaller than the maximum size permitted by the outer constraints, and we can fit its parameters.

Figure S1a shows re vs. for the same simulated fire data as that shown in Figure 5a-c. If the response was strictly exponential, we would find a straight line. Instead, we find some departure from exponentiality, so we can try to fit Equation (5). It follows from this equation that

(S1).

If we introduce Equation (5) into Equation (S1), we have

We define and , so

(S2).

Note that, in case of an exponential response, we would have =0, so  quantifies the difference from a simple exponential growth of as a function of re. Since the data in Figure S1a are distributed into groups corresponding to intervals of re, we can approximate the derivatives in Equation (S2) using differences between successive groups and thus:

(S3).

From Equation (S3), we expect a straight line if we represent vs. for our simulated data, as done in Figure S1b. I found a good fit to a straight line, except for the spot corresponding to the smallest , which was therefore ignored in the rest of the analysis (I discuss the behaviour of this spot below). I estimated  by simple regression, as shown in Figure S1b (the dashed line corresponds to  = 0).

Once we have , we can use it to extract the rest of the parameters with the help of the following equation, resulting from Equation (5):

(S4).

We must expect a straight line if we represent re vs. , as in Figure S1c. Once checked that this was the case (except for the spot with the smallest ), I estimated and k by simple regression. Since we know , we can extract k and . The set (k, , ) is easy to transform into the set (a, , ), so we have all of the parameters of Equation (5). The value rc can also be extracted from after calculating from the matrices of the simulated landscapes.

Once I obtained the parameters of Equation (5), I represented the expected response of to re in Figure S1a, along with the simulated response.

S4. Results for Section 5 (‘Analysis of the model using phase transition theory’)

The results agree with the hypothesis that the response of SOCFUS model simulated fire regimes to the environmental parameter re displays a second order phase transition. Figures S1b-c, corresponding to intermediate steps in the process of fitting the parameters, display straight lines as expected. After fitting the parameters, the response of the mean fire size to the environmental parameter re is well matched by Equation (5), in agreement with the phase transition hypothesis, as shown in Figure S1a. Only one of the spots, corresponding to the smallest re, is not fitted by this equation, and deviates from the straight lines in Figures S1b-c. This does not enable us to reject the hypothesis: the shape of Equation (5) is universally valid for second order phase transitions, but only as the critical threshold is asymptotically approached. Out of the ten spots in Figure S1, all are consistent with this asymptotic pattern except the one that is the furthest from the critical threshold. This result is completely consistent with our hypothesis.

The values fitted to the parameters of Equation (5) are: , , . This  is equivalent to . It differs slightly but unambiguously from  = 0 (dashed line in Figure S1b), which would correspond to an exponential response. Since 0.30, the above is equivalent to .

S5. Materials and Methods for Section 6 (‘Response of boreal forest fires to climate change: a forecasting exercise’)

S5.2.1. Weather index

I first tried to fit response functions to the results in Figure 5 and to use them for forecasting purposes. In principle, this would have enabled the p.d.f. to be predicted as a function of FWI. However, this led to a dead end, so I had to use SSR instead. I do not give the details of the methods that I used to fit the parameters of the response to FWI, which are quite extensive, and quite similar to those that I give for SSR below, but I give the reasons why I concluded that the use of FWI is inappropriate.

After fitting the response functions of  and sM to FWI for the period 1983-1995, I took the set of FWI measured for each of the fires of at least 1 ha in this period, calculated the p.d.f. expected for each of these FWI, and combined them into a single p.d.f. for each year. Then I compared the found and expected yearly  for 1 ha s < 512 ha, using non-parametric regression (Figure S2). The resulting slope was not one, as might have been expected, but 3.2. I carried out a nonparametric test of this result: (i) I calculated the residuals, i.e. the difference between found and expected  for each year; (ii) I calculated the squared Pearson’s correlation R2 between expected  and residuals; (iii) I randomised the residuals 9,999 times; (iv) I calculated R2 for each pseudoreplica; (v) I took both the original result and the results for each of the pseudoreplicas, and calculated the fraction of cases with R2 equal or larger than the original one, which is the significance level  that enables the null hypothesis of slope one to be rejected. I found = 0.015. Therefore, the slope is significantly different from one.

This discrepancy means that, when trying to make predictions for a given fire, the whole set of FWI in the year of the fire gives additional information not contained in the specific FWI registered in the closest station at the onset of the fire. This may be the result of correlations between FWI values in space and time, and of inertia in fuel moisture. Therefore, it is more appropriate to use a seasonal-scale weather index, such as SSR. There is an additional advantage to using this index: the SSR register extends beyond 1995, whereas FWI data cease to be available.

Alberta’s Forest Protection Division Weather Section kindly sent to me a series of the average SSR between 1983 and 2003 in a set of meteorological stations in the province.

For most of the analyses, I grouped the data in regular intervals of SSR. In this manner, all bins had an acceptable number of data. However, there were significant differences in this number, which needed to be dealt with appropriately (see below).

For the p.c.c.s. forecasts, I multiplied each of the values of SSR from 1983 to 2003 by 1.15, and used the resulting series as an input to the response functions that I found.

S5.2.2. Assumptions and their justification

S5.2.2.1. Shape of the p.d.f.

Following the results of Section 4, I assume that each possible state of fuel and atmosphere will be associated with a given p.d.f. of fire sizes, consisting of a bounded power law with specific parameters  and sM. If we assume that there is a simple cutoff at s = sM, this p.d.f. will obey Equation (3).

Since SSR is a measure that integrates the whole fire season, the p.d.f. that we find will integrate the p.d.f. corresponding to many different situations. As explained in Section 4.4, the combination of different power laws with similar  gives a p.d.f. that cannot itself be distinguished from a power law. However, the changes that take place during a year are often significant enough to produce a significant deviation from the power law. This deviation generally consists of a lower  at the upper range of the distribution, because high values of s mainly result from high FWI, which are associated with small . This effect is more apparent in years with highSSR: for the seasons with a low SSR, the statistical distribution of FWI is sharply peaked, while the variability in FWI increases with increasing SSR.

In principle, we could try to decompose one year’s p.d.f. into a linear combination of many power laws, in a manner analogous to a Fourier analysis. However, the deviation from a pure power law is too small and noisy to carry out such a decomposition. A much simpler approach seems to suffice for capturing all or almost all of the relevant information in our samples. This approach consists of separating the p.d.f. into two sections: a ‘lower’ section l, for s smaller than a given s, and an ‘upper’ section u, for s larger than s. There is a single  for each section, l and u, and a simple cutoff sM at the end of the upper section.