Why Small Fires Are Common and Big Fires Are Not Rare Enough

A Distribution With No Sensible Average

Out of the 120,000 fires in our main run, 13,373 of them burned exactly one cell and then went out. They make up 11.14% of every fire in the record. Between them they destroyed 13,373 cells out of 92,683,166. Rounded honestly, they destroyed nothing.

Now the other end: fires of a thousand cells or more made up 18.004% of the count and 88.5% of the burnt area, and the single largest of them took 38,758 cells, 14.8% of the whole lattice, once and never again. The median fire was 32 cells. The mean was 772.36.

A median of 32 against a mean of 772 states the whole problem. Whatever you want to call the typical fire here, the average is not it, and the typical fire is not the average either, because the tail is heavy enough that a handful of events almost nobody in the record ever lived through drags the arithmetic mean around behind them.

We started this because the model is charming. Two rules and nothing else, and out drops a size distribution shaped like the ones people plot for real wildfires and real earthquakes, which is exactly the charm and exactly the trap. We finished it less charmed, because the exponent everybody quotes for this model is not stable, our own goodness of fit test rejects the power law it is supposed to describe, and the reason the rejection does not show up more often in the literature is that a straight line drawn through a log-log plot cannot notice.

The Lattice Is Not a Forest

We have no forest. The club owns no plot, no drone and no fire permit, so our research is computation, and this article is the output of a program we wrote, running a mathematical object that is only distantly related to a wood. What we do have is a room with computers in it, and when we say we measured something, we mean we measured it from our own simulated fires, the way you would measure the frequency of heads by flipping a coin many times.

The model is due to Drossel and Schwabl in 1992 [2], building on the self-organised criticality programme that Bak, Tang and Wiesenfeld started five years earlier [1]. Each of the \(L \times L\) cells is empty or holds a tree, and we use helical boundary conditions: every cell carries a single index \(i\), its neighbours are \(i \pm 1\) and \(i \pm L\), and \(i + L^2\) is identified with \(i\). Grassberger wired his lattice the same way [7], so his densities and ours can be compared digit for digit.

Two moves, repeated forever:

• Growth. Pick a cell at random. If it is empty, a tree appears.
• Lightning. Pick a cell at random. If it holds a tree, the entire 4-connected cluster of trees containing it vanishes in one instant, and its size \(s\) is recorded as one fire.

Fires are instantaneous relative to growth, so nothing grows while anything burns. That double separation of time scales is what makes the model tractable, and it is what the review literature means by the critical state of this model [5]. It also means the two rates enter only through their ratio:

Each elementary event in our code is a growth attempt with probability \(\theta/(1+\theta)\) or a strike with probability \(1/(1+\theta)\), so the number of plantings between strikes is geometric with mean \(\theta\), while the critical limit is \(\theta \to \infty\), which is where the whole argument lives.

What the lattice leaves out is most of a forest, since fire spreads to every touching tree with probability one and to nothing else, so there is no wind, no slope, no ember thrown across a firebreak, no humidity and no species that resists burning. Trees appear at uniformly random cells, which no seed does, and burnt ground is ready for a new tree the instant the fire passes, so nothing in the model corresponds to soil damage, or to a seed bank, or to succession, or to the years a burnt slope spends as scrub. Lightning falls uniformly in space and time. Real ignition clusters in storms and follows drought. A cell is not an acre. A fire of 38,758 cells is not a hectare of anything, and we have no conversion to offer. Nothing in this article is a prediction about any real forest.

Checking the Ruler Before Measuring the Room

Before fitting anything from the model we checked the fitting machinery against data whose answer we already knew by construction. The exponents here sit near 1.15, well below 2, so the unbounded discrete power law is no use: \(\sum_{k} k^{-\tau}\) diverges for \(\tau \le 1\), leaving nothing to normalise with. We fit the power law truncated to the integer support \([s_{\min}, s_{\max}]\),

whose log-likelihood is \(\ell(\tau) = -\tau \sum_i \log s_i - n \log Z(\tau)\). Setting the derivative to zero gives a condition that is easy to state and easy to solve: the model's mean of \(\log s\) must equal the data's, and the model's mean decreases monotonically in \(\tau\), so one bisection finds the maximum. The standard error comes from the Fisher information, \(\mathrm{SE} = 1/\sqrt{n \operatorname{Var}(\log s)}\). The apparatus belongs to Clauset, Shalizi and Newman [12], with the truncation handled the way Deluca and Corral set out [13].

Then, synthetic data from a known exponent, handed to the estimator. Five tests, four of which used 200,000 draws:

Every \(|z|\) is below 2. The estimator recovers the truth inside its own stated error at two sample sizes, at three exponents and on two different supports, which is the whole test; had it failed, there would have been no article.

Then we fitted exactly the same synthetic data the way almost every undergraduate lab does it: histogram the values, take logarithms of both axes, fit a straight line by least squares. Figure 1 shows what that gives on data whose answer is known by construction.

The raw-histogram regression is not a little biased. At \(\tau = 1.5\) it returns 0.97215 when the answer is 1.5, an error of a third on data manufactured to have no complications at all. The mechanism is not mysterious: the tail bins hold one or two events each, empty bins are silently dropped because the logarithm of zero cannot be plotted, and the surviving points in the tail sit systematically above where the law says they should. The fit then flattens to accommodate them, and Newman set this out plainly twenty years ago, while regression through a log-log histogram is still the default method in a great many papers [12], including papers whose conclusions rest on the exponent it gets wrong.

Working Notes From the Club Table

The Arithmetic

One thing in this study is exact. Do it by hand before trusting anything that was fitted. In a steady state every tree that grows must eventually burn. Count per lightning strike. A strike is preceded by \(\theta\) growth attempts on average, of which a fraction \(1 - \rho\) land on empty cells and become trees, where \(\rho\) is the tree density. A fraction \(\rho\) of strikes hit a tree and start a fire. That fire has mean size \(\langle s \rangle\). Trees in must equal trees out:

No fitting, no approximation, no free parameter. Put our main run through it.

A fifth of a standard error. If the lattice were not in a steady state, or the flood fill were missing cells at the boundary, or the geometric gaps were off by one, these two numbers would part company and the whole study would be worthless. They do not part company, so the simulator is doing what the model says it should be doing. A check, not a discovery. Still, it is the only number here that comes with a guarantee.

The density trace backs it up. During the discarded warm-up of 10,485 fires the density averaged 0.39226. First half of the recording, 0.39187; second half, 0.39400. No trend worth naming survives 120,000 fires. A steady state is supposed to look exactly that boring.

The Test Says No

Now the distribution itself. Every fire in the main run, log-binned at eight bins per decade, with two lines drawn through it: the maximum likelihood fit over \(s \in [4, 500]\), and the line a least-squares regression puts through the raw log-log histogram over the same window.

It looks like a power law. The test says otherwise. We tested it with the semi-parametric procedure of Clauset, Shalizi and Newman [12]: fit the law, generate synthetic data sets of the same size from the fitted law, refit each one, and ask how often a synthetic set fits its own model worse than the data fit theirs.

Over the full range, \(s = 1\) to 38,758, the fitted exponent is \(\tau = 1.1174\), and the Kolmogorov-Smirnov distance between data and fit is 0.05978. Synthetic sets drawn from the fitted law produce a mean distance of 0.00166 and a 95th percentile of 0.00233. Our data sit 26 times the 95th percentile. Rejected, and not narrowly.

Over the fitted window \(s \in [4, 500]\) the fit is closer. Closer is not close enough: the distance of 0.00824 stands against a synthetic 95th percentile of 0.00356, and none of the 200 synthetic sets managed to fit its own model as badly as the data fit ours. \(p < 0.005\).

Before treating that as a result, look at what the same test does on subsamples of exactly the same fires, drawn from the same window and fitted the same way:

The fitted exponent barely moves; the verdict flips completely. The table says nothing about the forest fire model and everything about hypothesis tests: with enough data a test rejects any model that is even slightly wrong, and every model is slightly wrong. Somebody who ran 1,000 fires would publish that the power law fits beautifully and would be reporting the size of their sample, not a property of the world.

So the rejection on its own is cheap, and what makes it worth reporting is that we can see the shape of the failure rather than just its \(p\) value, and the shape is the one Grassberger reported as early as 1993, when high-statistics runs contradicted the exponents Drossel and Schwabl had derived [3], and which he pinned down in 2002 on lattices a thousand times larger than ours [7]. Pruessner and Jensen confirmed it with a better algorithm [8][10]. Figure 4 shows it. Multiply the distribution by \(s^{\tau}\). A true power law would come out horizontal.

The Exponent Will Not Sit Still

If the model were critical in the plain sense, \(\tau\) would be a number belonging to the model rather than to the run. We measured it at five values of \(\theta\) across two decades, with the window held at \([4, \theta]\) throughout so that it keeps the same position relative to the cutoff at every \(\theta\).

The exponent climbs from 1.0394 to 1.1741 across the table, and against the hypothesis that these five numbers estimate a single constant, the weighted mean comes to 1.1145 with \(\chi^2 = 466.3\) on 4 degrees of freedom, which is not a near miss. A constant exponent is not supported at any of the lattice sizes and ratios we can reach. Fitted as a straight line in \(\log_{10}\theta\), the drift is +0.069 per decade.

The errors in that table come from a moving-block bootstrap rather than from the likelihood, because consecutive fires in a single run are nothing like independent of one another, since a large fire clears the lattice and the fires immediately after it are small. Resampling contiguous blocks of the fire sequence keeps that structure. The main run barely notices: 0.0025 against the asymptotic 0.0027, a ratio of 0.92. Knowing that was worth the trouble.

One number in that figure is the best agreement in the whole study, and we want to be clear that it is also by some distance the cheapest. Our cutoff scaling exponent of 1.081 matches Grassberger's 1.08 to three digits, and it comes from a five-point regression whose own error bar is 0.019, which is to say from five points and a lot of goodwill. Agreement that good on that few points is partly luck, whereas the tree density comparison is the honest one: our \(\theta = 125\) run gives \(\rho = 0.38073\) from 262,653 lightning strikes, a binomial standard error of 0.00095, sitting 0.9 standard errors from Grassberger's tabulated \(0.379837 \pm 0.000006\) [7]. That comparison confirms, sharply and independently of anything we tuned, that we implemented the model he implemented.

Where the Lattice Runs Out

Two columns of that table should worry you. At \(\theta = 2000\) the largest single fire took 37.3% of the lattice. At \(\theta = 5000\) it took 59.0%, and a fire that eats three fifths of the world is not being limited by \(\theta\) any more, it is being limited by the edge of the simulation.

So we repeated the worst case on a lattice four times larger in area, and looked at what moved.

The exponent moved by −0.0078, which is 0.8 combined standard errors. Reassuring, and not proof: the cutoff proxy moved upward by 12.5%, which says the smaller lattice genuinely was truncating the tail, even though the exponent barely noticed that anything had gone missing. At \(\theta = 500\) the same check across \(L = 256\), 512 and 1024 gives 1.1175, 1.1269 and 1.1206, a spread of 0.009 against bootstrap errors of 0.005 to 0.006, so there is a wobble at the edge of significance and no trend.

Now the number that decides how much any of this is worth. Grassberger's 2002 study reached \(65536 \times 65536\) and \(\theta = 256{,}000\) [7]. Our largest lattice is \(1024 \times 1024\), and our largest \(\theta\) is 5,000, which is not close. We are short by 4,096 in area and by 51 in \(\theta\), and he concluded that the apparent scaling laws are transients, and that the true asymptotic regime might not arrive until \(\theta \sim 10^{40}\), a number with no physical home anywhere.

Take our own drift seriously for a moment and it says something faintly ridiculous: a straight line through our five points reaches \(\tau = 1.19\), Grassberger's estimate, at \(\theta \approx 10^{3.8}\), which is around 6,300, barely past where we stopped. If that extrapolation were right the argument would have ended in 2002, but the extrapolation is wrong, and the reason it is wrong is the reason to distrust extrapolations of this kind: the drift is decelerating, we cannot see where it stops, and a straight line fitted over two decades tells you nothing whatever about the twentieth decade. We print it only to show how easily a number like that fools you.

Plainly, then: our lattices are too small to settle whether this model is critical. We can confirm the drift. We can confirm that the compensated distribution is nowhere flat for long. We can confirm that a goodness of fit test rejects the power law, and we cannot tell a drift converging on 1.19 from one converging elsewhere, or from one that never converges, and no run we can afford would tell us. On the evidence we can generate, the honest position is the one Palmieri and Jensen describe [16]: the model is scale-free in appearance over the range anybody plots, and what it is doing asymptotically is a separate question that our hardware cannot touch.

The Strongest Objection to All of This

A reader who knows the field could put the whole article aside in one paragraph, and the paragraph would go roughly like this. Grassberger established in 2002, on lattices four thousand times larger, that this model does not show simple scaling, that the exponent drifts, and that the tree density had been badly underestimated. Pruessner and Jensen said the same with a different algorithm [8][10]. A club rerunning that on inadequate hardware has confirmed a known result and called it research. Nothing here is new.

We think that objection is mostly correct. What follows is our account of which parts survive it.

The drift itself is confirmation and we have labelled it as such: our \(\tau\) at \(\theta = 500\) sits 0.0231 below the pre-2002 value of 1.15 and our \(\tau\) at \(\theta = 5000\) sits 0.0241 above it, which is the drift restated in the form of a disagreement with a constant. Reproduction, not discovery.

The density agreement at \(\theta = 125\) is worth more than it looks, and Grassberger quotes \(\rho = 0.379837 \pm 0.000006\); we get 0.38073 with a standard error of 0.00095, putting us 0.9 standard errors away from him and nowhere near an accident. A number that sharp from an independent implementation, written from the published rules by people who had not seen his code, is the strongest evidence in this article that we are simulating the object he simulated. It also rules out the likeliest way for a student project to go wrong, which is to simulate a subtly different model and never notice the difference.

The part we would defend hardest is Figure 1 and the small table sitting beside it. The claim that log-log regression is a bad estimator of a power-law exponent is not new either [11][12][13], but the specific numbers are ours and they are checkable in a minute by anybody with a laptop: on 200,000 draws from a known \(\tau = 1.5\), maximum likelihood returns 1.50145 and least squares returns 0.97215. Most published fire-size exponents, we would guess, came out of the second method. Reed and McKelvey fitted North American fire records parametrically and found the power law inadequate over the full range [9]; Hantson, Pueyo and Chuvieco fitted global fire records by maximum likelihood and preferred a log-normal in most regions [14]; Corral and González, reanalysing several geoscience data sets with a careful protocol, found wildfires among the cases better described by a truncated log-normal [15]. Those three results are in tension with a large earlier literature that fitted straight lines, and Figure 1 is a small demonstration of why that tension exists.

A second objection we cannot dispose of. Our headline exponent depends on a window we chose. Across thirteen reasonable windows on the same 120,000 fires, \(\tau\) ranges from 1.0096 to 1.1390, a spread of 0.1294, while the statistical error on any one of them is between 0.0009 and 0.0049. Window choice moves the answer twenty-six times further than the largest of those errors, and fifty times further than the error on the value we report. We fixed the window before fitting and said so. The sensitivity is printed above, and that is the best response available to us, though none of it makes the underlying problem go away. When a distribution is not a power law, the exponent of the power law it is not becomes a property of your window rather than of the world.

What a Real Fire Season Would Have to Look Like

The toy stops being comfortable here, and here we have to be most careful, because this is exactly the place where a lattice of squares gets quietly promoted into a forest.

The interesting property of our distribution is not the exponent. What matters is that \(\tau\) sits below 2, because for \(P(s) \propto s^{-\tau}\) with \(\tau \le 2\), the mean of \(s\) does not converge as the upper limit grows; the integral \(\int s \cdot s^{-\tau}\,\mathrm{d}s\) diverges at the top. Our fires have a finite mean only because there is a cutoff, and the cutoff exists only because the lattice is finite and \(\theta\) is finite. Push the cutoff out and the mean goes with it, which is what the identity \(\langle s \rangle = \theta(1-\rho)/\rho\) says outright: the average fire is proportional to \(\theta\), so a world with less frequent ignition does not get smaller fires on average, it gets larger ones.

Read that sentence again with fire suppression in mind, and the arithmetic turns uncomfortable. In the model, reducing the ignition rate at fixed growth is exactly increasing \(\theta\). Our runs say what happens: mean fire size goes from 86.65 at \(\theta = 50\) to 7,433.18 at \(\theta = 5000\), a factor of 86 for a factor of 100 in \(\theta\), and the largest fire observed goes from 3,034 cells to 154,578. The total area burnt per unit time barely changes. In a steady state, everything that grows has to burn. What changes is how the burning is packaged: fewer, vastly larger events.

Now the cold water, and there is a lot of it. Our model has no suppression in it. Nobody fights these fires. It has no fuel age, no moisture, no wind and no climate, and the real mechanism by which suppression changes a fire regime runs through fuel accumulation and continuity, which is precisely what our uniform random planting fails to represent. Moritz and colleagues put the alternative on the table: real fire regimes, read off historical catalogues and a detailed fire spread model, look more like highly optimised tolerance than like self-organised criticality [11]. The difference bites on exactly the management question this paragraph keeps circling. Reed and McKelvey's fire records point the same way, and so do the global analysis of Hantson and colleagues and the careful reanalysis by Corral and González: real fire size distributions are frequently not power laws at all [9][14][15].

So what does our computation actually license anybody to say? Only this: one simple, fully specified system exists in which the average event size is set by the ratio of two rates, and in which suppressing ignition provably enlarges events. Its tail is fat enough that planning around the average is a category error, and whether any real region belongs to that class is a question about that region, and answering it needs fire records and fuel maps rather than a lattice. What the model can do is stop somebody being surprised. If your fire size distribution has an exponent below 2, the biggest fire you have seen is mostly a statement about how long you have been watching.

Where a Different Choice Would Have Changed the Answer

Several decisions in this study could have gone the other way, and some of them would have moved the headline number, and being specific about which is the point of this section.

The fitting window. That choice has the largest effect, by a distance. Fitting \([1, \theta]\) instead of \([4, \theta]\) gives 1.0526 rather than 1.1269. Fitting \([4, 5000]\), which runs past the cutoff, gives 1.0893. Fitting \([16, 500]\) gives 1.1390, and a reader who wanted the answer to be 1.15 could have got much closer by choosing \([6, 500]\), and a reader who wanted it to be 1.0 could have got there with \([64, 4000]\). Everything else hangs on one sentence: window fixed at \([4, \theta]\), chosen before fitting.

Maximum likelihood against regression. On our own fire data the two happened to agree to within 0.004, which is at once the luckiest result in the study and the most misleading one. On synthetic data the gap runs to 35%, as Figure 1 shows, and if we had used regression and got a plausible-looking number, we would have had no way to know whether it was the plausible number or the biased one.

Four neighbours against eight. Our fire spreads to the four orthogonal neighbours. Allowing the four diagonals as well raises the effective connectivity, which lowers the density at which large clusters appear and would change every density in this article. It would not change the qualitative picture, because the mean-field limit of the model is a critical branching process with \(\tau = 3/2\) in our convention [4] and the two-dimensional behaviour is well away from that, but every number in the tables would move.

Helical against periodic boundaries. Helical wrapping joins each row's right edge to the next row's left edge. The result is not quite a torus. We used it because Grassberger did, and we wanted the density comparison exact; on a \(512 \times 512\) lattice the difference is a single seam, and we would expect it to matter less than our bootstrap error, but we did not test that and so we cannot claim it.

Separated against unseparated time scales. The instantaneous-fire rule is the model, not an approximation to it, so relaxing it produces a different model in which fires and growth overlap and very large fires are partly extinguished by the gaps they create. That variant has different critical behaviour, and the published exponents do not describe it.

How much we ran. Our \(\theta = 5000\) run has 10,000 fires against 120,000 for the main run. Fires there cost about fourteen seconds per thousand. Its bootstrap error, 0.0056, is twice the main run's. More fires there would have sharpened the top of the drift curve, and would have done nothing at all about the fires hitting the lattice edge.

Reproducing this

One file, one command, Python 3.12 and numpy. No other package is used and nothing is downloaded.

python analysis/forest-fire-criticality.py > analysis/forest-fire-criticality-output.txt

Expected runtime is seven to twelve minutes, and the run printed in the output file took 699.3 seconds with numpy 2.4.2 on Python 3.12.3 under Windows 11, while a repeat on the same laptop with nothing else running took 425.7 seconds. The master seed is 20251124, hard-coded at the top of the file, and three independent streams are spawned from it through numpy's SeedSequence: one for the lattice, one for the fitting and bootstrap, one for the synthetic validation data. Peak memory is modest, a few tens of megabytes, because the lattice is a byte per cell and the largest temporary array is capped at 50,000 integers. The file pins BLAS to a single thread before importing numpy, which changes no number and lets the script start on a machine that is short of memory. We ran the finished file three times on the same machine and diffed the outputs: every scientific number matched to the last printed digit, and the only lines that differed were the wall-clock timings. If your output differs from ours by more than the printed standard errors, something is wrong somewhere and we would very much like to hear about it. The interactive model runs a smaller version of the same lattice in your browser, and at its default settings it reproduces the tree density and the fitted exponent to within 0.004, and the mass balance identity to a fifth of a percent.

References