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

A Distribution With No Sensible Average

Start at the small end. Out of 120,000 fires in our main run, 13,373 burned exactly one cell, which is 11.14% of every fire in the record, and between them those 13,373 fires destroyed 13,373 cells out of 92,683,166. Rounded honestly, they destroyed nothing.

Now the other end of the record, which is where the pleasure in this distribution actually lives. 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. The median fire was 32 cells; the mean was 772.36.

A median of 32 against a mean of 772. The whole problem of this distribution sits in the gap between those two numbers, and nowhere else. Whatever you want to call the typical fire in this model, the average fire is not it, and the average is not typical of anything, because the tail is fat enough to drag the arithmetic mean around by events that almost nothing in the record ever saw.

We started this because the model is beautiful. Two rules, nothing else, and out comes a distribution shaped like the ones people plot for real wildfires and real earthquakes; we finished the work admiring the model rather more and trusting the summary of it rather less. The exponent everybody quotes will not hold still, our own goodness of fit test rejects the power law it is supposed to describe, and the rejection stays out of the literature mostly because a straight line drawn through a log-log plot has no way of noticing.

The Lattice Is Not a Forest

The club owns no forest. No plot, no drone, no fire permit, and one room with computers in it. So our research is computation, and this article reports the output of a program we wrote, which runs a mathematical object only distantly related to a wood; 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 a great many times.

Drossel and Schwabl published the rules in 1992 [2], five years after Bak, Tang and Wiesenfeld had proposed that some systems tune themselves to a critical point without anybody turning a dial [1], and the forest became the friendliest example anyone had. The lattice is \(L \times L\) cells, each one either empty or holding a single tree. 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 our densities compare with his directly.

Criticality deserves spelling out once, patiently, since the rest of this article leans on the word. A physical system sits at a critical point when its correlation length has run away to the size of the system itself, so that a disturbance anywhere can reach anywhere, and the statistics of those disturbances stop having any characteristic size at all. Bak's claim was that certain systems arrive at that point on their own, with nobody holding a temperature steady to three decimal places, and the forest is the cleanest illustration anybody has drawn: growth pushes the lattice towards connectedness, fire knocks it back, and neither rule knows the other exists.

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 anywhere on the lattice while anything is burning. That double separation of time scales is the trick of the whole construction; it makes the model tractable, and it is what the review literature means when it calls this state critical [5]. The two rates therefore 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\). The critical limit sits at \(\theta \to \infty\), and every argument in this field lives there, which is unfortunate, because no computer ever built has been anywhere near it.

What the lattice leaves out is, more or less, everything that makes a forest a forest. 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 is behaviour no seed from a real tree has ever shown. 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 arrives in storms, follows drought, and cares a great deal where it is. A cell is not an acre, a fire of 38,758 cells is not a hectare of anything, and nothing in this article is a prediction about any real forest.

Checking the Ruler Before Measuring the Room

Nothing here was fitted before the fitting machinery had been tested on data whose answer we already knew by construction. The exponents in this model sit near 1.15, well below 2, so the unbounded discrete power law is no use to us at all: its normalising sum \(\sum_{k} k^{-\tau}\) diverges for \(\tau \le 1\), and a distribution that cannot be normalised cannot be fitted. We labour this because the habit of fitting first and validating never is what produced a literature of exponents nobody can reconcile, and because the check costs an afternoon. 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 as Deluca and Corral set out [13], and none of it is ours; what is ours is the decision to test it before using it.

Then synthetic data, drawn from an exponent we picked ourselves and kept out of the estimator's reach. Five tests, four of them with 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: bin the values into a histogram, take logarithms of both axes, then fit a straight line by least squares. Here is what that method gives on data whose answer we know by construction rather than by argument.

The raw-histogram regression is not a little biased; it is wrong in a way that matters. At \(\tau = 1.5\) it returns 0.97215 when the answer is 1.5, an error of a third, on data with no noise in it beyond the sampling. 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. Newman set this out plainly twenty years ago, and the method is still the default in a great many papers [12], which is the part that annoys us.

Working Notes From the Club Table

The Arithmetic That Has to Balance

One quantity in this study is exact. Do it by hand, before trusting anything that came out of a fit. 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 land on a tree. Each of those starts a fire of mean size \(\langle s \rangle\). Trees in must equal trees out:

No fitting and no free parameter anywhere in it. Put our main run through it.

A fifth of a standard error. If the lattice were not in a steady state, if the flood fill were missing cells at the boundary, if the geometric gaps were off by one, or if the warm-up were too short, these two numbers would part company and the whole study would be worthless. They do not, so the simulator is doing what the model says, which is a check rather than a discovery, and no other number in this article arrives with a guarantee attached to it.

The density trace backs it up: during the discarded warm-up of 10,485 fires the density averaged 0.39226, over the first half of the recording it averaged 0.39187 and over the second half 0.39400, and no trend worth naming shows up across 120,000 fires, which is exactly what a steady state is supposed to look like.

The Test Says No

Here is the distribution. 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 on the plot, and it is not one under the test. 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 is 0.05978, against a synthetic mean distance of 0.00166 and a synthetic 95th percentile of 0.00233. Our data sit 26 times the 95th percentile, rejected and not narrowly. The distribution that looked so convincing on a log-log plot does not survive being compared against data drawn from its own fitted law, which is the only fair comparison there is.

Inside the fitted window \(s \in [4, 500]\) the agreement is far better and the verdict does not change: the distance is 0.00824 against a synthetic 95th percentile of 0.00356, and not one of 200 synthetic sets did worse than the data. \(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 flip 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. 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], and Figure 4 shows it. Multiply the distribution by \(s^{\tau}\). A true power law of that exponent would flatten into a horizontal line. It does not flatten. Every curve leans, and then lifts into a shoulder just below its cutoff, and the shoulder is the interesting part of this whole study, because a pure power law has no room in it for a feature with a size.

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. Against the hypothesis that these five numbers estimate one constant, the weighted mean is 1.1145 and \(\chi^2 = 466.3\) on 4 degrees of freedom, so a constant exponent is not supported at any size we can reach. An exponent that depends on \(\theta\) is not an exponent belonging to the model; it is a reading taken at a particular distance from a limit nobody can reach, and the pleasure of that result is that it turns a memorised constant back into an open question. 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 not independent of one another. A large fire clears the lattice and the fires immediately after it are small. Resampling contiguous blocks of the fire sequence keeps that structure. In the main run it makes almost no difference, 0.0025 against the asymptotic 0.0027, a ratio of 0.92, and knowing that the correlations cost us so little was itself worth the trouble.

One number in that figure is the best agreement in the whole study, and we would like to be clear that it is also the cheapest one in it. The cutoff scaling exponent 1.081 against Grassberger's 1.08 matches to three digits, out of a five-point regression whose own error bar is 0.019, and agreement that good on that few points is partly luck. 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]. Sharp, independent confirmation that we implemented the model he implemented. Getting that number right is the least glamorous thing in the article and the one we were happiest about, because it means the object on our laptop is the object in the literature rather than a cousin of it.

Where the Lattice Runs Out

The last two columns of the \(\theta\) table are the ones that should make a reader uneasy, and they made us uneasy for a fortnight. At \(\theta = 2000\) the largest single fire took 37.3% of the lattice. At \(\theta = 5000\) it took 59.0%. 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.

The exponent moved by −0.0078, which is 0.8 combined standard errors. Reassuring, and not proof. The cutoff proxy moved by 12.5%, upward, which tells you that the smaller lattice genuinely was truncating the tail even though the exponent barely noticed the truncation. 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 ran lattices up to \(65536 \times 65536\) and reached \(\theta = 256{,}000\) [7], while our largest lattice is \(1024 \times 1024\) and our largest \(\theta\) is 5,000. We are short by a factor of 4,096 in area. In \(\theta\) we are short by a factor of 51. His conclusion was that the apparent scaling laws are transients, and that the true asymptotic regime might not arrive until \(\theta \sim 10^{40}\), a number we enjoy looking at. A regime that begins at \(\theta \sim 10^{40}\) is not a regime anybody will ever simulate, which is either a devastating objection to the whole self-organised criticality story or a statement about how patient physics has to be, and we cannot tell you which.

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. The extrapolation is wrong, and the reason it is wrong is the reason to distrust every extrapolation 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 the number only to show how easily one fools you.

So, plainly. Our lattices cannot settle it. We can confirm that the exponent drifts, that the compensated distribution is nowhere flat for long, and that a goodness of fit test rejects the power law. We cannot distinguish a drift that converges to 1.19 from one that converges somewhere else, and we cannot rule out one that never converges at all. 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 a small version of that, finding a drift, and announcing it has confirmed a known result on hardware too small for the question. Nothing here is new. The objection is one we put to ourselves before anybody else could, and the honest reply is not a defence of the whole article but a sorting of it into the part that merely reproduces and the part we would still put our names to if Grassberger had never published a line.

The objection is mostly correct. We want to say 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. A reproduction, not a discovery.

The density agreement at \(\theta = 125\) is worth more than it looks. Grassberger quotes \(\rho = 0.379837 \pm 0.000006\); we get 0.38073 with a standard error of 0.00095, which puts us 0.9 standard errors from a number measured on lattices we could never run. 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 find out.

The part we would defend hardest is Figure 1 and the table 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. We would guess the second method is behind most reported fire-size exponents. 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. The window choice moves the answer by twenty-six times the largest of those errors, and by fifty times the error on the value we actually report. We fixed the window before fitting and reported the sensitivity, which is the best available response, and it does not make the underlying problem go away. When a distribution is not a power law, the exponent of the power law it is not turns out to be a property of your window rather than of the world.

What a Real Fire Season Would Have to Look Like

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

The exponent itself is not the interesting property of this distribution, whatever value it takes. What matters is that \(\tau\) sits below 2. Below 2 is where a distribution stops having a typical member, and the whole of the discomfort in this section comes from that one inequality rather than from anything about trees. 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 out with it, obediently, and without any limit at all. The identity \(\langle s \rangle = \theta(1-\rho)/\rho\) says exactly this: 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 it is 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, because in a steady state everything that grows must burn; what changes is how the burning is packaged, into fewer and vastly larger events.

Now the cold water, of which there is a great deal more than we would like. 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 argued, using both historical catalogues and a detailed fire spread model, that real fire regimes are better described by highly optimised tolerance than by self-organised criticality, and that the difference matters for exactly the management question this paragraph is circling [11]. 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. There exists a simple, fully specified system in which the average event size is set by the ratio of two rates and in which suppressing ignition provably enlarges events. The tail is fat enough that planning around the average is a category error. 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 honestly do is stop somebody being surprised by a fire they had no name for. 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, some of them would have moved the headline number, most would not, and being specific about which is more useful than a general promise of care. None of what follows is a confession; it is the list a reader needs in order to work out how much of our headline number is a fact about the Drossel-Schwabl model and how much is a fact about us.

The fitting window. Nothing else in the whole study moves the headline answer anything like as far as this. 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. 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 in this article hangs off a single sentence in the working notes, namely that the window was fixed at \([4, \theta]\) before we had seen one fitted number.

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 in it. On synthetic data with a known answer they disagree by up to 35% (Figure 1). 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 of a burning cell and to nothing else. 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 the right edge of one row to the left edge of the next, which makes the lattice a long spiral rather than a torus. We used it because Grassberger did, and the density comparison had to be 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. The 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, because a thousand fires there cost about fourteen seconds of wall clock. 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. How much of that you get depends entirely on what else the machine is doing. 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 directly in your own browser. 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