The emergence of scale-free fires in Australia

The emergence of scale-free fires in Australia arXiv:2110.10014v3 [cond-mat.stat-mech] 16 Feb 2023 Giorgio Nicoletti,1, 2 Leonardo Saravia,3, 4 Fernando Momo,5, 6 Amos Maritan,1 and Samir Suweis1 1 Laboratory of Interdisciplinary Physics, Department of Physics and Astronomy “G. Galilei”, University of Padova, Padova, Italy 2 Department of Mathematics “T. Levi-Civita”, University of Padova, Padova, Italy 3 Centro Austral de Investigaciones Cientı́ficas (CADIC - CONICET), Ushuaia, Argentina 4 Área Biologı́a y Bioinformática, Universidad Nacional de General Sarmiento, Buenos Aires, Argentina 5 Universidad Nacional de General Sarmiento, Instituto de Ciencias, Área Biologı́a y Bioinformática, Buenos Aires, Argentina 6 Instituto de Ecologı́a y Desarrollo Sustentable (INEDES - CONICET), Buenos Aires, Argentina Between 2019 and 2020, during the country’s hottest and driest year on record, Australia experienced a dramatic bushfire season, with catastrophic ecological and environmental consequences. Several studies highlighted how such abrupt changes in fire regimes may have been in large part a consequence of climate change and other anthropogenic transformations. Here, we analyze the monthly evolution of the burned area in Australia from 2000 to 2020, obtained via satellite imaging through the MODIS platform. We find that the 2019-2020 peak is associated with signatures typically found near critical points. We introduce a modeling framework based on forest-fire models to study the properties of these emergent fire outbreaks, showing that the behavior observed during the 2019-2020 fire season matches the one of a percolation transition, where system-size outbreaks appear. Our model also highlights the existence of an absorbing phase transition that might be eventually crossed, after which the vegetation cannot recover. INTRODUCTION Bushfires are an intrinsic part of Australia’s landscape dynamics. Its natural ecosystems have evolved to coexist with fires, and mitigation strategies to reduce their impact have been learned in the most vulnerable areas [1]. Yet, the 2019-2020 fire season was particularly catastrophic. It began in July 2019, at the end of the country’s hottest and driest year on record, and wildfires were unprecedented in their spatial extent and severity [2–5]. In the eastern Australia states of New South Wales and Victoria, around 5.8 million hectares of mainly temperate broadleaf forest were burned by a series of high-impact fires, many of which exceeded a size of 100, 000ha and continued to burn for weeks after ignition. Several studies highlighted how this abrupt departure from the historical trend may have been in large part a consequence of climate change and other anthropogenic transformations [5–10]. Furthermore, these high-impact fires had a devastating effect on Australia’s biodiversity. Of more than 830 taxa - comprising birds, reptiles, frogs, mammals, and freshwater fish - around one-fourth lost to the fires between 10% and 50% of their Australian habitat, sixteen of them lost between 50% and 80%, and three more than 80% [4]. These drastic changes, with their catastrophic effects on vegetation and on biodiversity, are often associated with critical transitions, i.e., conditions that inevitably lead to large-scale fire outbreaks and subsequent widespread damage [11–14]. Such behavior has been observed in many different systems ranging from Amazon forests [15, 16] to Kalahari vegetation [17] and more in general in tropical forests fragmentation [18, 19]. In physical systems with many degrees of freedom, these phenomena are well-known to appear at the edge of phase transitions. When a system undergoes a continuous phase transition at a critical point, scale-free behaviors described by power-laws are found - such as long-range correlations and diverging susceptibility to external perturbations - due to the underlying scale-invariance that emerges at criticality [20–23]. This lack of a characteristic scale is a possible mechanism behind the abrupt appearance of large and out-of-scale events, such as the high-impact fires experienced by Australia between 2019 and 2020. In this work, we analyze the monthly evolution of the burned area in the East and Southeast temperate broadleaf and mixed forests of continental Australia [24]. These data, spanning from November 2000 to June 2020, are obtained via satellite imaging through the MODIS platform [25], and allow us to analyze the spatiotemporal properties of fire propagation. Unsurprisingly, we find that the 2019-2020 peak of the burned area exceeds the historical data. Then, thanks to the high spatial resolution of the data, we study the distributions of spatially-separated clusters of burned area, as well as their evolution in time. By applying tools from Statistical Physics, we find that during 2019-2020 the distribution of fire outbreak sizes is compatible with a power-law, and it is invariant under spatial coarse-graining. Our results suggest that such fires lacked a characteristic size, and thus the system may have been poised at a critical point of their spreading dynamics [22]. To understand the drivers and the type of such critical transitions, we introduce a paradigmatic spatial model that describes the concurrent spreading of fires and vegetation over a two-dimensional lattice. In a regime where the timescale of fire propagation is much faster than the vegetation one, our numerical simulations suggest that the model 2 FIG. 1. The cumulative distribution of the fire sizes at different years. (A) The timeseries of the normalized number of burned pixels per month in the East and Southeast temperate broadleaf and mixed forests of continental Australia, Nburned (t), from 2000 to 2020, normalized by the total number of pixels. Years are defined as the twelve months occurring between June and May. The year 2019-2020 largely exceeds the peaks of the previous twenty years. (B-E) For a given year, we can compute the cumulative fire size distribution on a nearest-neighbor basis. Even though peaks, such as during 2002-2003, often display either longer tails in the distribution or are dominated by few, very large fires, a distinctive power-law behavior emerges during 2019-2020. See also Figure S2. predicts the crossing of a percolation-like transition [26] to a more arid climate, where spreading becomes easier for the fires and harder for the vegetation. Differently from self-organized forest-fires models [27–33], the dynamics of our model depends only on two effective ecological parameters. When these parameters cross the percolation-like critical point, the features of the model - such as the distribution of the fires’ sizes - are qualitatively comparable to the ones observed during the bushfire season of 2019-2020 in Australia. This suggests that this kind of phase transition in the vegetation-fires dynamics may have been at the heart of the emergence of scale-free fire outbreaks. Our paradigmatic model encompasses another kind of critical point as well, that corresponds to an absorbing phase transition [34, 35] after which the vegetation cannot recover. Although possibly unrelated from an ecological perspective, it foreshadows how critical points may lead to further abrupt and fundamental changes in the fire-vegetation dynamics. RESULTS Fire sizes distributions In order to shed a light on how the 2019-2020 bushfire season emerged, we analyze the timeseries of the burned area obtained from 236 monthly satellite images of Australia, spanning from November 2000 to June 2020 (Methods and Figure S1). Each month is represented by a binary matrix Mt , whose entries (Mt )ij represent an area of 500m2 and are set to 1 if the corresponding pixel matches an area that has burned in the span of that month. The exceptional nature of the 2019-2020 events Pis perhaps already striking from the timeseries of the total burned area spanning the last 20 years, Nburned (t) = ij (Mt )ij , plotted in Figure 1a. As a gauge of the extent of the damage, less than half of the pixels were burning during the second largest peak - which took place in the season 2002-2003. Most importantly, given the spatial nature of our data, we can also compute the cumulative distribution of clusters of burned pixels for a given month - on a nearest-neighbors basis - starting from Mt . Hence, for each month (i) nc (t) and each matrix Mt , we compute the number of clusters nc (t) and their sizes {Ac }i=1 . We identify nc with the (i) number of separate fire outbreaks, and Ac with the corresponding outbreak sizes. In particular, we can compute the corresponding cumulative fire sizes distribution during a given year (Methods), defined from June to May to include 3 FIG. 2. Properties of the timeseries of the number of fires and their maximum size. (A) Plot of Nfires (t), the number of fire outbreaks in a given month, and of Mfires (t), the largest outbreak of a given month. Both timeseries have been normalized by their maximum value in order to compare them. (B) The phase of the Hilbert transform of both the timeseries of the number of fires and of the maximum fire size show a major synchronization during 2019-2020. (C) Similarly, the corresponding amplitudes suggest that during the 2019-2020 fire season a very large number of fires coexisted together with extremely large ones. This behavior is indeed captured by the power-law behavior of the fire size distribution. (D) The Kuramoto index between the phases of the Hilbert transform of nc (t) and of mc (t) and the correlation between the respective moduli. A clear synchronization emerges during 2019-2020, and in the same year the correlation between the moduli spikes as well. This behavior is compatible with the power-law cumulative distribution of the fire sizes found in the same year. the summer of the Southern Hemisphere. These distributions, shown in Figure 1b-e, typically display longer tails during peaks of the burned area Nburned (t), whereas the sizes are exponentially suppressed if the overall burned area is low (see also Figure S2). In particular, higher peaks of Nburned (t) - e.g., the 2002-2003 or the 2019-2020 season - are associated with distributions that span a wide range of sizes. Although it is tempting to relate such distributions to power-laws, the finite size of the system and the limited data make detecting such power-laws a non-trivial task [36–39]. Hence, it is paramount to understand whether the distribution of 2019-2020, characterized by a cutoff that is typically associated with the finite size of the system (Figure 1c), is quantitatively different from the ones of previous years. Notably, the range of fire sizes of previous years - e.g., 2002-2003 or 2006-2007 - show that fires larger than the ones in 2019-2020 took place, see Figure 1b and Figure 1d. Therefore, to gain further insights into the fire dynamics, we extract the timeseries of the number of fires per month Nfires (t) = nc (t) - i.e., the number of connected clusters of each matrix Mt - and of the size of the largest fire Mfires (t) (i.e., the size of the largest connected cluster of each matrix Mt ). In general, we do not expect these timeseries to be synchronized and, indeed, we typically see that a high number of clusters in a given year usually does not imply large clusters as well, as we can see in Figure 2a. To study the relation in time between Nfires (t) and Mfires (t), we perform a dynamical analysis of these two timeseries by introducing the phase and the modulus of their Hilbert transform (Methods). We plot them in Figure 2b-c. Then, for every year, we compute the associated Kuramoto index [40], which measures the synchronization between the number of fires and the maximum fire sizes, and the moduli correlation (see Figure 2d and Methods for further details). We find that during the 2019-2020 season Australia experienced neither the largest fire in our data nor the largest number of fire outbreaks - but rather a major synchronization between the two timeseries emerges. That is, both the number of fires and the size of the largest fires suddenly increased with respect to previous years. Such presence of many and very large outbreaks at once can be interpreted as a distinctive proxy of the widespread damage the outbreaks caused. Furthermore, it suggests that the distribution we see in 2019-2020 may be associated to features that we expect to see in power-law distributions and scale-free phenomena. Crucially, such a scale-free distribution might be tightly related to the dramatic impact that the 2019-2020 bushfire season had on the vegetation and on biodiversity. In fact, fire sizes that are distributed as a power-law give rise to both a few large fire clusters - corresponding to the distribution’s long tails - and many 4 smaller ones - the bulk of the power-law distribution. As the former devastated entire regions, the latter created pockets of vegetation fuel in other areas that could possibly act as an ignition to the next high-impact fire. Such a catastrophic departure from the historical trend suggests that a fundamental shift in the underlying dynamics might have occurred. In the next section, we will test the hypothesis that the 2019-2020 distribution is compatible with a scale-free one by introducing a spatial coarse-graining. Spatial coarse-graining The spatial resolution of our data allows us to carefully test the hypothesis that the power-law distribution we see during 2019-2020 is a signature of an underlying scale-invariance. Here, we draw inspiration from the Renormalization Group concepts [20, 41–43] and implement the so-called coarse-graining (CG) step to understand if such scale-invariance is present. We perform, at each time, a spatial CG through a block-spin transformation of Mt , by grouping together nearby pixels in 2 × 2 plaquettes. We then define the new super-pixels through a majority rule, in such a way that if the plaquette contained a majority of burned pixels, the corresponding coarse-grained pixel will be burned as well, and vice-versa (Methods). Then, we follow the properties of the system along these repeated transformations. In fact, a coarse-graining transformation amounts to studying a system at different spatial scales. If the system is truly scale invariant, we expect that its properties will not change under repeated CG steps. Hence, and compatibly with the quality of the data, if the distribution of the fire size is a true power-law it will remain a power law after one or more CG transformations, with a corresponding finite-size scaling correction (Methods). In principle, one should iterate the coarse-graining indefinitely, to unravel the properties of its fixed points - however, with real data we are limited by the finite size of our system. Since each of the coarse-graining steps we are employing reduces the linear size of the system by half, after four CG steps we are left with a matrix that contains only ≈ 0.4% of the initial number of pixels. If only few but large fires are present in the original system, this coarse-grained version will be dominated by system-sized outbreaks. On the other hand, if many but small fires characterized the initial state, the coarse-graining transformations will drive the system to a configuration where virtually no fires are present. In particular, the behavior of probability distributions along the coarse-graining is particularly relevant in determining the properties of a critical system [44]. In Figure 3a-c we show an example of the effect of three coarse-graining steps on a sub-region of the matrix P M = t Mt where each entry indicates the total number of times the corresponding pixel has been burned. As we can see, the spatial coarse-graining preserves some features of the original matrix, although the number of pixels is reduced by a factor 24 . As a reference, in Figure 3d we show once more the timeseries of the burned area Nburned and in Figure 3e-h we follow the cumulative distributions of the fire sizes along the coarse-graining for selected years. For instance, during the 2007-2008 season or the 2015-2016 season the CG quickly suppresses the fire sizes distribution and only small fires are left. Interestingly, if we consider a season associated with a marked peak of the burned area, such as the 2002-2003 season in Figure 3e, we see that, although the distribution keeps its distinctive long tails along the coarse-graining, the bulk of the distribution changes and no evident cut-off appears. On the other hand, and crucially, during the 2019-2020 season the bulk of the power-law distribution of the fire sizes is left invariant, while the cut-off associated with successive CG steps is poised at smaller and smaller system’ sizes (see Methods). Quantitatively, in Figure 3i we show that maximum-likelihood fits [36] of the distributions at different levels of coarse-graining in 2019-2020 the exponents remain compatible with one another. Then, we repeat the same procedure for the other years and plot the corresponding exponents in Figure 3j. Taking into account the standard deviation of these exponents, only the distribution associated with the 2019-2020 season consistently displays a power-law behavior with the same exponent at all CG steps. These results strongly support the power-law nature of the 2019-2020 distribution, revealing a unique underlying scale-invariance in the spatial structure of the fire outbreaks taking place in that season. This scale-invariance, in turn, manifests itself dynamically as a synchronization between the number of outbreaks and their maximum size. Altogether, our data suggest that the emergent properties we observe are related to a phase transition. Thus, a fundamental question arises: what has driven the 2019-2020 fire dynamics close to what appears to be a critical point? Paradigmatic model for the vegetation-fires dynamics To qualitatively understand the abrupt changes observed during the 2019-2020 bushfire season, we introduce a minimal stochastic model of the forest-fire class [28]. Differently from classical forest-fire models that display self-organized criticality [31, 45, 46] and multistability [47], we describe the concurrent stochastic spread of both fires and vegetation 5 FIG. 3. The properties of the data under spatial coarse-graining. (A-C) An example of the effect of three coarsegraining steps on the overall number of fires per pixel in a sub-region of our data. The more steps are performed, the less the number of pixels left. (D) The timeseries of the burned area Nburned plotted in Figure 1, shown as a reference. (E-H) The coarse-graining corroborates the presence of a very robust scale invariance during 2019-2020, whereas in the previous years the shape of the distribution is significantly changed by the CG transformation. For instance, during 2007-2008 and during 2015-2016 the distribution of the fire sizes is exponentially suppressed, and after four coarse-graining steps there are almost no more fires to begin with. The distribution of 2002-2003, although not exponential, lacks the expected finite-size cutoffs and displays changes in the bulk of the distribution along the coarse-graining. (I) Power-law fit of the coarse-grained distributions in 2019-2020 using maximum-likelihood fitting methods. (J) If we fit a power-law to each year, we find that 2019-2020 displays the most consistent exponents at different coarse-graining steps. For comparison, in, e.g., 2006-2007 the exponents vary in the range ≈ [1.5, 1.8], which is almost four times larger than the one found in 2019-2020. See also Figure S3. between neighboring nodes of a network. Although extensions of the forest-fire model with different vegetation growth and with climate effects have been proposed [48, 49], our approach seeks to include only minimal features to understand whether they are able to explain the patterns observed in our data. Without fires, the vegetation V is free to spread to its nearest neighbors at a rate λV on a given graph - for instance, a 2-dimensional lattice - and spontaneously disappears with a death rate dV . Then, a fire F can ignite on a vegetation site with rate bF and spreads with a rate λF over an effective topology that is determined by the structure of the vegetation clusters. At the same time, the 6 FIG. 4. The fundamental properties of our model. (A-B) A depiction of the model dynamics as a multi-layer graph and the corresponding transition rates. (C) On a 2D lattice the model displays a charge-discharge behavior if the vegetation dynamics is much slower than the fire one, and fires are relatively rare events. Here (dF , bF , λF ) = (25, 10−5 , 500) and (dV , λV ) = (0.5, 3). (D-E) The vegetation layer undergoes an isotropic percolation transition at ζVperc ≈ 2.63 where a spanning cluster appears. In (D) we plot the size of the largest vegetation cluster cmax and in (E) the mean vegetation cluster size χV , which peaks V at the transition. Both the plots are from a 250 × 250 lattice (see Figure S4). (F-G) If we consider fires that spread over a fixed vegetation configuration (Methods), (F) below the percolation threshold ζV < ζVperc the cumulative distribution of the fire size sF is always exponentially suppressed due to the small vegetation clusters. Above it (G), the fires may spread on a spanning cluster, and therefore we have system-size outbreaks if ζF is small enough. (H-I) Comparison of the analytic solution of the mean field equations and a stochastic simulation on a fully connected network with 500 nodes. (H) With parameters (dF , bF , λF ) = (1, 0.5, 10) and (dV , λV ) = (0.1, 0.5) the absorbing state, i.e., the empty configuration is the stable mean-field solution. (I) For (dF , bF , λF ) = (10, 0.1, 100) and (dV , λV ) = (1, 3), instead, noise-induced oscillations around the mean-field stationary values emerge. Notably, the mean-field approximation is not able to predict the charge-discharge behavior described above, which is a consequence of spatial effects that are thus fundamental in our model. vegetation cannot occupy a site with a fire F , thus both the topology of the fire layer and of the vegetation layer change dynamically with time. Once a fire is over, with a rate dF , the corresponding site will become an empty site ∅ for the vegetation layer, and will not be present in the fire layer. Hence, although archetypal, our model is described in terms of few parameters that can be thought of as functions of environmental conditions. Notably, a similar model was proposed by Zinck et al. [50] to analyze data from the Canadian Boreal Plains. However, the modeling approach proposed by the authors did not include nearest-neighbor spreading nor a death rate for the vegetation. As we will see, in our model, these parameters are fundamental in shaping the spatial structure of fires. Indeed, heuristically, we can think of our model as defined on a multi-layer network [51–53]. In this depiction, 7 the topology of the vegetation layer is fixed, but the vegetation sites dynamically govern the topology of the fire layer, as we sketch in Figure 4a. Hence, we expect the interplay between the spatial spreading of both vegetation and fires to be a crucial feature of our model, whose rates are shown in Figure 4b. The vegetation alone obeys λ V Vi + ∅j∈∂i −−→ Vi + Vj d (1) V Vi −−→ ∅i , where i is a site and ∂i is the set of the neighbors of i. These reactions for the vegetation dynamics, thus correspond to the well-known contact process [34, 35, 54, 55], an archetypal model of absorbing phase transitions. We highlight that, differently from most SOC and previous models, we do not include an immigration term for the vegetation i.e., an external field in the contact process. This amounts to assuming that vegetation can only spread from other vegetation sites, rather than reappearing in random sites. On top of this dynamics the fire spreading is determined by the reactions λ F Fi + Vj∈∂i −−→ F i + Fj b F Vi −−→ Fi (2) d F Fi −−→ ∅i . These reactions, if considered independently from Eq. (1), represent instead a contact process with resource depletion - meaning that the empty sites are unavailable for fires to spread. Model simulations and time-scale separation We perform exact stochastic simulations of the model on a 2-dimensional lattice using the Gillespie algorithm (Methods) [56]. Crucially, for the model to be reasonable, we must assume that the vegetation dynamics is much slower than the one of the fires and that the birth rate of the fires bF is typically very small [50]. In Figure 4c we show that the model in this range of parameters indeed displays a charge-discharge dynamics, with long periods of almost undisturbed vegetation spreading followed by shorter periods of fire spreading following the rare ignition of an outbreak. This time-scale separation limit corresponds to the assumption that the vegetation configuration does not change during the propagation of a fire. Therefore, we study how a fire propagates on top of a fixed stationary vegetation configuration. In this scenario, the phase space is described by the adimensional parameters ζF = dF /λF and ζV = λV /dV (Methods). A small value of ζF gives rise to fires that are extremely effective at spreading and, viceversa, a large value of ζV implies a quick vegetation regrowth. Remarkably, since the vegetation layer in the absence of fires follows a simple contact process, we expect a percolation transition at ζVperc ≈ 2.6, as recently shown with numerical simulations [57]. At this value, a system-size cluster of vegetation appears, coexisting with a significant number of distinct, but smaller vegetation clusters (Figure 4d-e). At ζV = ζVperc , in particular, if L → ∞ an infinite cluster of vegetation appears. If we call cV the vegetation cluster size and n(cV ) the number of vegetation clusters of size cV , we can define the mean vegetation cluster size ratio between the first two moments as P 2 cV (cV ) n(cV ) χV = P (3) cV cV n(cV ) where the sum runs over all vegetation clusters. This quantity, which we plot in Figure 4e, is expected to diverge at the percolation transition due to the scale-free In fact, in percolation theory, P Pnature of the vegetation cluster sizes. this is nothing but the mean cluster size χ = cV cV wcV , with wcV = cV n(cV )/ cV cV n(cV ) the probability that a vegetation site belongs to a cluster of size cV , which displays a power-law divergence close to the percolation threshold [26]. This transition has a crucial impact on the cumulative distribution of the fire sizes sF , as we see in Figure 4f-g. In fact, below the percolation transition of the vegetation, fires are severely limited by the size of the vegetation clusters, and thus the distribution of sF is exponentially suppressed even at small ζF . On the other hand, above the percolation threshold, the vegetation clusters tend to be larger, and fires can be large if ζF is small enough. This suggests, as highlighted in [50], that a critical transition may underlie the vegetation-fires dynamics. Although this percolation-like transition emerged from the spatial nature of our model, we can also solve it analytically in a mean-field approximation - which amounts to ignoring such spatial features to begin with (Methods). Yet, the mean-field solution allows us to reveal the presence of yet another critical point, an absorbing phase transition 8 FIG. 5. The properties of the time-scale separated model and its behavior under spatial coarse-graining. (A) At a given set of parameters (ζV , ζF ) we plot the ratio between the mean fire size hsF i and the mean size of a vegetation cluster hcV i in a 250 × 250 lattice. The black dotted lines represent contour lines of χ̃ (the product of the number of vegetation clusters ncV and the maximum fire size smax ), which is maximized around the percolation transition ζVperc for low enough values of F 5 0 ζF . (B-G) We seed nF = 10 fires on a lattice with linear size L = 1000 in order to study the distribution of the fire sizes sF and the corresponding coarse-grained distributions. (B-C) At low values of ζV , the cumulative distribution of the fire sizes is exponential and is further suppressed along the coarse-graining at all values of ζF . (D-E) At ζVperc , if ζF is low enough the fire size distribution becomes a power-law that is invariant along the coarse-graining. (F-G) For high values of ζV , on the other hand, the system is dominated by few large clusters of vegetation, and the corresponding large fires are highlighted by the coarse-graining. This regime is not particularly realistic at low ζF , since it would require climate conditions that allow for large fires, i.e., a warm and arid climate, but at the same time for an extremely effective vegetation spread. [34, 35, 55]. This phase transition separates a phase in which the mean-field stable configuration predicts a non-zero density of both fire and vegetation from a phase in which the stable configuration is the empty one, see Figure 4h-i. Crucially, the mean-field picture is drastically different from a spatially embedded model. Indeed, the spatial structure significantly changes the way fires spread due to the modulation of the underlying vegetation structure, leading to isotropic percolation which will play a fundamental role. Model coarse-graining and the emergence of scale-free fires Simulations of the model allow us to study the properties of the area burned by fires at different values of (ζV , ζF ). In particular, in Figure 5a we show the behavior of the ratio between the average fire size hsF i and the average vegetation cluster size hcV i observed in configurations with given parameters (ζV , ζF ). This parameter is fundamental because it helps us understand the potentially damaging effects of the fires on the underlying vegetation substrate. Whenever hsF i / hcV i ≈ 1, it implies that a fire that originates in a given vegetation cluster has a non-vanishing probability to burn the entire cluster. In Figure 5a we plot another relevant quantity as well - the black dotted lines represent the contour lines of χ̃ = smax × ncV , where ncV is the number of vegetation clusters. This quantity is particularly significant because ncV F can be interpreted as a rough estimate of the number of possible fires in the system, whereas smax tells us how large F they can be. In the data, these two quantities both reached high values at the same time during 2019-2020. 9 In order to study the behavior of the fire sizes distribution under the same spatial coarse-graining applied in the data, we choose n0F = 105 fire seeds in a large lattice of linear size L = 1000 (Methods). Then, we analyze the resulting burned area in a given point of the phase space (ζV , ζF ). In particular, we look at the distribution of the fire sizes and, thanks to the large size of the lattice, at how it changes along repeated CG transformations. We find four different regimes, shown in Figure 5a-g. If ζV is high enough, typically the vegetation can spread effectively and regrow any burned vegetation. Yet, if ζF is low, fires can propagate almost unboundedly due to the underlying large vegetation clusters. The resulting distributions, shown in Figure 5g), are therefore dominated by very large fires. Indeed, since for ζV > ζVperc a spanning cluster is present, vegetation sites that are far away are likely connected and fire can spread from one to the other. This regime is perhaps unrealistic since it leads to extremely large fires in an otherwise vegetation-rich environment. Yet, a similar dynamics was observed in fire-prone communities where species with post-fire recruitment have the most flammable canopies [58]. On the other hand, a more realistic regime is described by a high ζV and a high ζF as well. This regime corresponds to environmental conditions that favor a vegetation-rich system and suppress fires, and therefore we expect to see a small burned area. In fact, as we see in Figure 5f, fires are small as they are not able to propagate effectively, not even on the underlying spanning cluster of vegetation sites. Crucially, in both these regimes (Figure 5f-g) the coarse-graining accentuates the tails of the fire size distribution, since the coarse-graining will unravel the largest fires that propagate on the vegetation spanning cluster. On the other hand, if ζV is low, vegetation regrowth is typically suppressed. In this case, when ζF is high, fires tend to be small as we see in Figure 5b, but so do the clusters of vegetation. Indeed, hsF i / hcV i can dangerously increase because substantial parts of the underlying vegetation clusters can burn even at high ζF . Finally, when ζF is also low, not only the vegetation clusters can hardly regrow, but a fire can systematically burn the entire cluster in which it originates since hsF i / hcV i ≈ 1. This regime is not sustainable in the long time - the fires are likely to outpace the vegetation regrowth and eventually desertification will take place. Notably, this regime cannot be distinguished from the distribution of the fire sizes alone, in Figure 5c, which stays exponential due to the lack of large vegetation clusters. The vegetation percolation transition lies in between these regimes, and it is here that power-law distributed fires emerge at low enough ζF (Figure 5d-e). In fact, at this point we see a distinctive scale-invariant configuration emerging following a power-law distribution with an invariant bulk under spatial coarse-graining. Moreover, this is also the region where max sF × ncV is maximized, since the system can experience large fires coexisting with a large number of clusters of vegetation. Therefore, the features that we observe during the 2019-2020 fire season are best described by our model in a time-scale separation approximation close to the vegetation percolation transition ζV = ζVperc , with small fire suppression ζF . Note that in this regime, hsF i / hcV i remains small. However, the mean itself is not representative in the presence of scale-free distributions, hence it is not a reliable index of fire damage anymore. DISCUSSION How did Australia reach such a critical point? Although our modeling approach is paradigmatic, it provides a clear physical interpretation of its control parameters ζV and ζF . Indeed, their value is determined by climate conditions. Therefore, prolonged droughts, higher temperatures, and a more arid climate - all recognized as contributors to the 2019-2020 bushfire season [6, 9, 10] - might have pushed both ζV and ζF to lower and lower values, eventually reaching and crossing the percolation transition between 2019 and 2020. Notably, the 2019-2020 year has been unusually hot and dry in part due to natural meteorological phenomena, such as a shift in the polar winds above Antarctica and one of the strongest positive swings in the Indian Ocean Dipole. The former contributed to stratospheric warming, which in turn contributed to bringing hot, dry weather to much of Australia. The latter, in its positive phase, may have led to a reduction in rainfall over the southern and most northerly regions of Australia [6]. However, on top of - and possibly as a cause of - all this natural variation, global warming is making the country even hotter and drier [59], with the devastating effects that we highlighted in this work. Finally, the mean-field analysis and the vegetation layer of our paradigmatic model predict the presence of yet another critical point, one of a very different nature associated with the absorbing phase transition of the contact process [34, 35, 54, 55] at λV /dV := ζVabs ≈ 1.6. This phase transition separates a phase in which the only stable configuration is the absence of vegetation, and a phase in which vegetation is present [60]. Crucially, with the addition of the fire dynamics, a slow enough vegetation spreading implies that fires at high values of ζF can burn large clusters of vegetation. This scenario may push the system to a state in which the vegetation goes extinct. Such states are much harder to reach in more realistic and highly complex dynamics of fire spreading in forests. For example, one should consider that broadleaf Australian forest species, such as Eucalyptus, have resilience and resistance traits, like re-sprouting and seed banks, that allow for a rapid post-fire recovering even in intense fire-regimes [61, 62]. Yet, repeated fires with short return times would cause the exhaustion of these capacities [60]. These considerations 10 do suggest that the isotropic percolation transition observed during the 2019-2020 bushfire crisis may foreshadow a worsening condition that, in the far future, might push the system to a forest-savanna-like type of transition [63, 64]. Overall, our results suggest that the unprecedented bushfire season that Australia experienced between 2019 and 2020, with outbreaks appearing at all scales, is compatible with a phase transition in the vegetation-fires dynamics driven by a worsening climate. Our work shows how phase transitions and critical points play a fundamental role in shaping this dynamics, and their presence and consequences will be more and more relevant as climate change will quickly deteriorate the climatic conditions. LIMITATIONS OF THE STUDY Future works should aim to develop quantitative methods to infer the values of the model’s parameters from data, both from both fires spreading and vegetation evolution. Although the present study lacks such inference steps, procedures such as simulation-based inference [65] may be well-suited to this aim. In particular, ecological and environmental drivers evolve over time, both due to seasonality and climate change. This would amount to prescribe a dynamics for the parameters ζF and ζV of our model, as well as bF , which are instead considered constant in our analysis. Changes in these parameters over time may affect the dynamics, creating feedback effects that are taken into account in the present study. Notably, quantities from Information Theory may be a promising extension to disentangle the environmental effects from the vegetation-fire interaction [66]. It will be crucial to account for and disentangle contributions coming from natural variations and from the anthropogenic impact, in order to assess mitigation strategies that are becoming more and more vital. Furthermore, here we only apply from the Renormalization Group in the form of coarsegraining and finite-size scaling. It will be of particular interest to consider other phenomenological approaches to the Renormalization Group, beyond simple coarse-graining [67]. Finally, it will be paramount to apply the analysis carried out in this work to other areas of the world where large and extended fire outbreaks are appearing. ACKNOWLEDGMENTS AM was supported by “Excellence Project 2018” of the Cariparo foundation. LS was financially supported by ”Fondo para la Investigación Cientı́fica y Tecnológica (FONCYT) PICT 2020 SERIE A 2628”. AUTHOR CONTRIBUTIONS L.S. and S.S. designed the study. L.S., A.M., and S.S. supervised the research. L.S. and F.M. collected the data. G.N. analyzed the data. G.N., A.M., and S.S. designed and studied the model. G.N. performed the simulations and the analytical calculations. G.N. prepared the figures. G.N. and S.S. interpreted the results and wrote the first draft of the manuscript. All authors contributed to the final version of the article and approved the submitted version. DECLARATION OF INTERESTS The authors declare no competing interests. DATA AND CODE AVAILABILITY This paper analyzes existing, publicly available data from the NASA MODIS platform. All original code has been deposited at Zenodo and is publicly available at [68, 69]. 11 METHODS Data collection We defined the region of study as the East and Southeast temperate broadleaf and mixed forests of continental Australia using the ecoregions defined by Dinerstein [24], accessible at http://ecoregions2017.appspot.com/, which represents an area of 48 · 106 ha (see Figure S1). For this region, we estimate the burned areas using the NASA Moderate-Resolution Imaging Spectroradiometer (MODIS) burnt area Collection 6 product MCD64A1 [25], which is a monthly product with a 500m pixel resolution. We downloaded the images, using Google Earth Engine, as geoTIFF and then we converted them to a binary matrix (circa 4000x8000) using the R statistical language [70]. Then, for each month we have a binary matrix Mt , whose pixels represent an area of 500 m2 and can be either 1 - if there has been a fire in that pixel in the span of that month - or 0 - if no event occurred, meaning that no burned area was detected. Cluster distributions We define a cluster of a binary matrix Mt using a nearest-neighbors connectivity, i.e. the pixels that belong to a cluster are defined using the connectivity matrix   0 1 0 (4) Cbasis = 1 1 1 0 1 0 which defines the usual nearest neighbors of a 2-dimensional lattice. We also repeated the analysis described in the main text using a next-nearest-neighbors connectivity and the results do not change significantly. Therefore, for each (i) nc (t) . Then the cumulative matrix Mt we end up with a number of clusters nc (t) and the areas of each cluster {Ac }i=1 fire size distribution of Mt can simply be evaluated as   nc (t) θ A(i) − s X c F P (sF ) := P (Ac > sF ) = (5) nc (t) i=1 where θ(·) is the Heaviside function. To evaluate yearly distributions, we pooled the cluster sizes from all matrices Mt of a given year. This amounts to assuming that the clusters found in subsequent months are independent. Indeed, we find that the overlap between burned pixels in Mt and Mt+1 is always small, with respect to the number of burned pixels. Cluster dynamics We can exploit the timeseries of both the number of clusters and their areas to probe the underlying properties of the fire dynamics. In particular, we look at the number of clusters nc (t) = Nfires (t) and the area of the largest cluster (i) nc (t) mc (t) = maxi {Ac }i=1 = Mfires (t). We normalize both these timeseries by dividing them by their maximum value, in order to make them comparable (Figure 2a). Other normalizations, such as a standard z-score, give essentially the same results. In order to understand how the evolution of these two timeseries relates in time, we introduce the Hilbert transform of a real-valued timeseries x(t) as Z ∞ i x(t + τ ) − x(t − τ ) dτ (6) H[x(t)] = x(t) + lim π ε→0 ε τ which is a complex timeseries. Thus we can compute its phase ϕx (t) = arctan Im[H[x(t)]] Re[H[x(t)]] and its modulus ρx (t) = q Im2 [H[x(t)]] + Re2 [H[x(t)]] and how they change in time (Figure 2b-c). We can further quantify the relations between nc (t) and mc (t) by looking at the Kuramoto index [40] of their Hilbert transforms and at the correlation between the corresponding moduli. We define the Kuramoto index on a given year as D E Kyear = eϕnc (t)−ϕmc (t) (7) year 12 and the correlation between the moduli as hρnc ρmc iyear − hρnc iyear hρmc iyear Cyear = r h i. Q 2 2i hρ − hρ i i year i year i∈{mc ,nc } (8) In Figure 2d we do indeed see that Kyear becomes significantly close to 1 during 2019-2020, hence the two timeseries are highly synchronized during the year. Similarly, the correlation between the moduli has a positive spike in the same period. It is worth noting that this is true even if in the original timeseries neither the number of clusters nor the size of the largest one are maximal during 2019-2020. The fundamental change in the behavior of the system is that, during this year, both of them peak in a synchronized fashion, which leads to the power-law distribution shown in the main text. Spatial coarse-graining A quantitative and powerful way to assess the scale-invariance of a system is given by a properly defined coarsegraining procedure [20, 41, 43]. In the spirit of Statistical Physics, a suitable coarse-graining for a binary matrix Mt is a block-spin transformation of the associated 2-dimensional square lattice. Namely, the k-th coarse-graining step (k+1) (k) amounts to define a super-pixel σi′ from the previous pixels σi via the majority rule ( P (k) > ⌊card(Bi )/2⌋ 1 if (k+1) j∈Bi σi (9) σ i′ = 0 otherwise where ⌊·⌋ is the floor function and Bi is the i-th set of pixels such that {Bi } forms a non-overlapping covering the original 2-dimensional lattice. In particular, we take Bi ∈ M(2 × 2) so that at each coarse-graining step the number of pixels is reduced to a fourth of the original ones and therefore we can perform enough coarse-graining steps. Notice P (k) that, in this case, the majority rule is not exact since the cardinality of Bi is even. Thus, if j∈Bi si = 2 we (k+1) randomly assign the value of si′ to be either 0 or 1. In the spirit of the Renormalization Group, we should follow physical observables and - in particular - probability distributions [44] to look for scale-invariance along the coarse-graining. That is, if the system is scale-invariant in a spatial sense we should see that, even if we are coarse-graining the system, some of its properties will not change up to some finite-size cutoff, because the small-scale features are indistinguishable from the large-scale ones. This is exactly what we look for when we compare the cumulative probability distributions of the cluster sizes at different coarse-graining steps. As at each coarse-graining step we observe a smaller and smaller system, we can exploit finite-size scaling. Thinking of a percolation-like transition [26], the probability distribution of the fire sizes in a system of linear size L scales as s  F −τ +1 Pcumulative (sF ) = sF ψ (10) LD where D is related to the critical exponent of the correlation length and τ is the exponent of the power-law distributed fire sizes. In particular, D is the fractal dimension of the fires. Hence, for a properly chosen value of D, we expect that Pcumulative (sF )sτF−1 as a function of sF /LD will collapse onto the same curve. We find this collapse with D ≈ 1.95, which suggests once more, and in terms of the Renormalization Group, that the 2019-2020 fire seasons appear to behave like a system close to a phase transition (see Figure S3). In fact, the fractal dimension tells us the size sF of a fire outbreak changes with its linear size, i.e., sF ∼ LD . Notice that in bond percolation we would expect the fractal dimension to be D = 91/48 ≈ 1.896 in a two-dimensional lattice, which is compatible with what we find in the data. However, the exponent τ of the fire size distribution is different from the one expected in bond percolation, suggesting that the universality class might be different. Let us note that, in our model, bond percolation only happens in the isolated vegetation layer, and not in the layer where fire propagates. Contact process and critical points The contact process [34, 35, 54, 55] is an archetypal model for absorbing phase transitions, which describes spreading phenomena over a set of sites {σi }i=1,...,N . Each site can be either occupied σi = 1 or empty σi = 0. Empty sites are 13 occupied by neighboring occupied sites at a rate λ, whereas occupied sites become empty at a rate µ. The mean-field equations for the density of occupied sites ρ is given by ρ̇ = ρ(λ − µ) − λρ2 . This equation has two stationary solutions. The first one is the empty configuration ρvst = 0, which is only stable if λ < µ. The empty configuration is an absorbing configuration, that is, once it is reached the system cannot leave it since no reactions are possible. If λ > µ, the stable stationary solution is ρast = 1 − µ/λ. The value λabs := µ is the critical point of the model at which the absorbing phase transition takes place - below λabs , the system will always reach the absorbing empty configuration, whereas above λabs non-zero values of ρ are possible. Conversely, ρ is the order parameter of the system, which identifies the two phases. This kind of critical point is present in our model of vegetation-fire spreading as well, as we can show analytically in the mean-field case. The contact process on the 2D lattice, however, displays another kind of phase transition related to its spatial structure [57], a percolation transition. The order parameter of this transition is the probability that a site belongs to a spanning cluster, i.e., an infinite cluster, which is zero below the percolation transition and greater than zero above. Notice that Martı́n and collaborators [57] use a slightly different definition of the contact process, in which empty sites are occupied by neighbors with a probability p and occupied sites become empty with a probability 1 − p. They show numerically that the percolation transition in a 2D lattice happens at pperc ≈ 0.725. We can immediately recover our formulation by noting that p = λ/(λ + µ), giving the result used in the main text (λ/µ)perc ≈ 2.63. Mean-field behavior of the model The mean-field equations read dp∅ = −λV pV p∅ + dF pF + dV pV dt dpF = −dF pF + (λF pF + bF )pV dt dpV = −(dV + bF + λF pF )pV + λV pV p∅ dt (11) where p∅ , pF and pV are the probabilities of each state. Since p∅ = 1 − pF − pV we consider only the equations for pF and pV . In general, the stationary state of the system is given by ( dF pF = (λF pF + bF )pV . (12) (dV + bF + λF pF )pV = λV pV (1 − pV − pF ) abs These equations have an absorbing solution, since (pabs V , pF ) = (0, 0) is a trivial solution of the system. Importantly, it is easy to show that by adding a birth term for the vegetation this empty solution disappears, as expected. The abs Jacobian matrix evaluated at (pabs V , pF ) is given by   −dF bF (13) J abs = 0 −(bF + dV ) + λV = λV −bF −dV . Thus, the empty state is only stable below λabs = −dF and µabs whose eigenvalues are µabs 2 1 V = bF +dV , which is the absorbing critical point of the system. Notice that λF does not play a meaningful role in the stability of the empty state, a feature that is likely wrong in a spatially embedded model. The other stationary state of the system is given by q stat −dV λF − (bF + dF + λF )λV + fF,V stat pF = (14) 2λF λV q stat λF (−2bF − dV ) − (bF + dF − λF )λV + fF,V (15) pstat = V 2λF (λF + λV ) stat 2 abs st where fF,V = 4dF λV λF (λabs V − λV ) + (dV λF − (bF + dF + λF )λV ) is positive above λV . The eigenvalues of J always have a negative real part if λV > bF + dV while they always may have a non-vanishing imaginary part. Hence, the relaxation towards the steady state typically happens in an oscillatory fashion. In particular, these oscillations play a major role in the evolution of the finite-size stochastic model, where the noise can push the system to the absorbing state or produce sustained stochastic oscillations, as we see in Figure 4i. 14 Exact stochastic simulation Simulations of the three-state model on a given network, such as a 2-dimensional lattice, are performed using the Gillespie algorithm [56]. If we assume that there are N sites in the network and M possible transitions - in our model, (t) M = 6 - then, at each time the network can be associated with a propensity matrix Aµi , where µ = 1, . . . , N and (t) i = 1, . . . , M . Each row of Aµi is given by the transition rates that the µ-th site can undergo, given its state at time P P (t) (t) t. We introduce the total propensity α0 = µ i Aµi , so that the waiting time for the next transition is given by (t)
−1 τ (t) = − α0 log u (16) where u is uniformly distributed in [0, 1]. Then, the transition ī that occurs and the site µ̄ at which it occurs are such that µ̄−1 ī−1 XX (t) (t) Aµi ≤ α0 v < µ̄ X ī X (t) Aµi (17) µ=1 i=1 µ=1 i=1 where v is once again uniformly distributed in [0, 1]. We then update Aµ̄i with the new transition rates for µ̄ and set the time to t + τ . Simulations and time-scale separation As pointed out in the main text, it is reasonable to expect that the vegetation dynamics is much slower than the fire dynamics. However, the parameter space of the model is extremely large and thus a phase space plot for the full model proves to be unfeasible. Therefore, in order to simplify the problem and reduce the number of free parameters, we assume that the vegetation configuration does not change during the propagation of a fire. This approximation is compatible with the time evolution of the model if we also assume that fires are rare events and the complete model predicts a charge-discharge dynamics. In fact, realistically, we expect the vegetation - if dry enough - to act as fuel during a fire propagation, which thus has to stop when locally all fuel is exhausted. Then the vegetation regrows and only after enough fuel is accumulated a new fire can start - that is, the two processes happen at different
timescales. This assumption is implemented in simulations as follows. Let us start with a network {σi }N i=1 , {Eij }) where the N sites σi are such that σi ∈ {∅, V } and {Eij } are the edges between the sites. We look for a stationary configuration {σistat }N i=1 of the reactions λV ,Eij (σi = V ) + (σj = ∅) −−−−→ (σi = V ) + (σj = V ) (18) dV (19) (σi = V ) −−→ (σi = ∅) λV ,Eij where the notation −−−−→ means that the reaction happens at a rate λV if and only if i and j are joined by an edge Eij . The system has an absorbing configuration {σi = ∅} and its stationary configurations only depend on the ratio of the reaction rates ζV = λV /dV . See Supplementary Figures for examples of such configurations in a 2-dimensional lattice. Our approximation consists in obtaining a network over which the fires can spread from the stationary configuration {σistat }N i=1 . In particular, we consider the subgraph induced by the map g : i 7→ µ defined for all the indexes  i such V that σi = V . If we call these sites sµ = σg(i) , we end up with the vegetation subgraph {sµ }N µ=1 , {Eµν }) where Eµν = Eg(i)g(j) and NV is the number of original vegetation sites. This subgraph is typically composed of many disjointed components. These components contain roughly the same number of nodes for ζVabs < ζV ≪ ζVperc , since the stationary configuration is dominated by a large number of small vegetation clusters, whereas as we approach ζVperc a giant component emerges and it is eventually dominant for ζV ≫ ζVperc . We now assume that sµ ∈ {∅, V, F }, and notice that the initial configuration is such that sµ = V , ∀µ = 1, . . . , NV . In order to sample the distribution of the fire sizes we can choose a site sµ̄ and set sµ̄ = F . Then, to simulate a fire, we consider the reactions λF ,Eµν (sµ = F ) + (sν = V ) −−−−−→ (sµ = F ) + (sν = F ) d F (sµ = F ) −−→ (sµ = ∅) (20) (21) 15 until there are no more F sites in the network. Thus, the fire dynamics only depend on ζF = dF /λF . The fire size i.e., the burned area - is simply the number of empty sites N∅ of the final configuration. One should be careful that if sµ̄ is chosen at random between all sites we typically favor larger components of the vegetation subgraph. Thus, we first uniformly sample a given component Cs of the vegetation subgraph, and then we randomly choose a site within the selected component and set Cs ∋ sµ̄ = F . This assumption is qualitatively equivalent to the assumption that if two fires start in the same cluster they will contribute to the same burned area. To be precise, this is only true if ζF is small enough, so that two fires inside the same component will coalesce with high probability. However, for larger values of ζF we expect fires to be small independently of the size of the underlying component, hence our assumption does not affect the results. In this way, we are now able to computationally explore the model’s behavior effectively and systematically. In particular, for each value of ζV , we simulate a large number of stationary configurations {σistat }N i=1 . Then, for each of these configurations, at a given value ζF we simulate a number of fires much larger than the number of components Cs , thus ending up with a set of burned areas {N∅ } that gives us the fire size distribution at (ζF , ζV ). 1 SUPPLEMENTARY FIGURES FIG. S1. Region of study. The shaded area represents the region of study encompassing the East and Southeast temperate broadleaf and mixed forests of continental Australia. FIG. S2. Comparison of yearly distributions of fire sizes from 2001-2002 to 2019-2002. 2 FIG. S3. Invariance under coarse-graining and data collapse. (a) The cumulative distribution of the fire sizes during 2019-2020 along the coarse-graining transformations, with the best power-law fit obtained by maximum-likelihood [36]. The exponents remain compatible at different CG steps. Notice that the different τ reported here those of the fire size distribution, not of the cumulative distribution. (b) The cumulative distributions of the data at different CG steps collapse into the same curve, once appropriately rescaled, as predicted by finite-size scaling. In order to appreciate the quality of the collapse, notice the different units in the vertical axis with respect to panel (a). (c) The two parameters τ and D that are needed to collapse the cumulative distributions of the fire sizes can be chosen so that the collapse error will be minimal [71]. As expected, we find τopt ≈ 1.65, which is compatible with the fitted exponent. FIG. S4. Stationary configurations of the vegetation contact process. Different colors represent clusters of different sizes in a 250 × 250 2-dimensional lattice and t different values of ζV . Notice how at the percolation threshold a system-size cluster appears. 3 [1] C. P. Yates, A. C. Edwards, and J. Russell-Smith, “Big fires and their ecological impacts in australian savannas: size and frequency matters,” International Journal of Wildland Fire, vol. 17, no. 6, pp. 768–781, 2008. [2] D. B. Lindenmayer and C. Taylor, “New spatial analyses of australian wildfires highlight the need for new fire, resource, and conservation policies,” Proceedings of the National Academy of Sciences, vol. 117, no. 22, pp. 12481–12485, 2020. [3] P. Deb, H. Moradkhani, P. Abbaszadeh, A. S. Kiem, J. Engström, D. Keellings, and A. Sharma, “Causes of the widespread 2019–2020 australian bushfire season,” Earth’s Future, vol. 8, no. 11, p. e2020EF001671, 2020. [4] M. Ward, A. I. Tulloch, J. Q. Radford, B. A. Williams, A. E. Reside, S. L. Macdonald, H. J. Mayfield, M. Maron, H. P. Possingham, S. J. Vine, et al., “Impact of 2019–2020 mega-fires on australian fauna habitat,” Nature Ecology & Evolution, vol. 4, no. 10, pp. 1321–1326, 2020. [5] N. Phillips and B. Nogrady, “The race to decipher how climate change influenced australia’s record fires,” Nature, vol. 577, no. 7791, pp. 610–613, 2020. [6] A. J. Dowdy, “Climatological variability of fire weather in australia,” Journal of Applied Meteorology and Climatology, vol. 57, no. 2, pp. 221–234, 2018. [7] P. Yu, R. Xu, M. J. Abramson, S. Li, and Y. Guo, “Bushfires in australia: a serious health emergency under climate change,” The Lancet Planetary Health, vol. 4, no. 1, pp. e7–e8, 2020. [8] M. M. Boer, V. R. de Dios, and R. A. Bradstock, “Unprecedented burn area of australian mega forest fires,” Nature Climate Change, vol. 10, no. 3, pp. 171–172, 2020. [9] N. B. Arriagada, D. M. Bowman, A. J. Palmer, and F. H. Johnston, “Climate change, wildfires, heatwaves and health impacts in australia,” in Extreme Weather Events and Human Health, pp. 99–116, Springer, 2020. [10] N. J. Abram, B. J. Henley, A. Sen Gupta, T. J. Lippmann, H. Clarke, A. J. Dowdy, J. J. Sharples, R. H. Nolan, T. Zhang, M. J. Wooster, et al., “Connections of climate change and variability to large and extreme forest fires in southeast australia,” Communications Earth & Environment, vol. 2, no. 1, pp. 1–17, 2021. [11] M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara, “Early-warning signals for critical transitions,” Nature, vol. 461, no. 7260, pp. 53–59, 2009. [12] E. H. van Nes, B. M. Arani, A. Staal, B. van der Bolt, B. M. Flores, S. Bathiany, and M. Scheffer, “What do you mean,‘tipping point’ ?,” Trends in ecology & evolution, vol. 31, no. 12, pp. 902–904, 2016. [13] J. C. Rocha, G. Peterson, Ö. Bodin, and S. Levin, “Cascading regime shifts within and across scales,” Science, vol. 362, no. 6421, pp. 1379–1383, 2018. [14] M. Scheffer, Critical transitions in nature and society. Princeton University Press, 2020. [15] B. Wuyts, A. R. Champneys, and J. I. House, “Amazonian forest-savanna bistability and human impact,” Nature Communications, vol. 8, no. 1, pp. 1–12, 2017. [16] T. E. Lovejoy and C. Nobre, “Amazon tipping point,” Science Advances, vol. 4, no. 2, p. eaat2340, 2018. [17] T. M. Scanlon, K. K. Caylor, S. A. Levin, and I. Rodriguez-Iturbe, “Positive feedbacks promote power-law clustering of kalahari vegetation,” Nature, vol. 449, no. 7159, pp. 209–212, 2007. [18] F. Taubert, R. Fischer, J. Groeneveld, S. Lehmann, M. S. Müller, E. Rödig, T. Wiegand, and A. Huth, “Global patterns of tropical forest fragmentation,” Nature, vol. 554, no. 7693, pp. 519–522, 2018. [19] L. A. Saravia, S. Doyle, and B. Bond-Lamberty, “Power laws and critical fragmentation in global forests,” Scientific Reports, vol. 8, p. 17766, 2018. [20] J. J. Binney, N. J. Dowrick, A. J. Fisher, and M. E. Newman, The theory of critical phenomena: an introduction to the renormalization group. Oxford University Press, 1992. [21] G. Caldarelli, R. Frondoni, A. Gabrielli, M. Montuori, R. Retzlaff, and C. Ricotta, “Percolation in real wildfires,” Europhysics Letters, vol. 56, pp. 510–516, nov 2001. [22] D. Sornette, Critical phenomena in natural sciences: chaos, fractals, selforganization and disorder: concepts and tools. Springer Science & Business Media, 2006. [23] D. McKenzie and M. C. Kennedy, “Power laws reveal phase transitions in landscape controls of fire regimes,” Nature Communications, vol. 3, no. 1, pp. 1–6, 2012. [24] E. Dinerstein, D. Olson, A. Joshi, C. Vynne, N. D. Burgess, E. Wikramanayake, N. Hahn, S. Palminteri, P. Hedao, R. Noss, M. Hansen, H. Locke, E. C. Ellis, B. Jones, C. V. Barber, R. Hayes, C. Kormos, V. Martin, E. Crist, W. Sechrest, L. Price, J. E. M. Baillie, D. Weeden, K. Suckling, C. Davis, N. Sizer, R. Moore, D. Thau, T. Birch, P. Potapov, S. Turubanova, A. Tyukavina, N. de Souza, L. Pintea, J. C. Brito, O. A. Llewellyn, A. G. Miller, A. Patzelt, S. A. Ghazanfar, J. Timberlake, H. Klöser, Y. Shennan-Farpón, R. Kindt, J.-P. B. Lillesø, P. van Breugel, L. Graudal, M. Voge, K. F. Al-Shammari, and M. Saleem, “An Ecoregion-Based Approach to Protecting Half the Terrestrial Realm,” BioScience, vol. 67, pp. 534–545, June 2017. [25] L. Giglio, W. Schroeder, and C. O. Justice, “The collection 6 MODIS active fire detection algorithm and fire products,” Remote Sensing of Environment, vol. 178, pp. 31–41, 2016. [26] D. Stauffer and A. Aharony, Introduction to percolation theory. CRC press, 2018. [27] P. Grassberger and H. Kantz, “On a forest fire model with supposed self-organized criticality,” Journal of Statistical Physics, vol. 63, no. 3, pp. 685–700, 1991. [28] B. Drossel and F. Schwabl, “Self-organized critical forest-fire model,” Physical review letters, vol. 69, no. 11, p. 1629, 1992. [29] P. Grassberger, “On a self-organized critical forest-fire model,” Journal of Physics A: Mathematical and General, vol. 26, no. 9, p. 2081, 1993. 4 [30] K. Christensen, H. Flyvbjerg, and Z. Olami, “Self-organized critical forest-fire model: Mean-field theory and simulation results in 1 to 6 dimenisons,” Phys. Rev. Lett., vol. 71, pp. 2737–2740, Oct 1993. [31] B. D. Malamud, G. Morein, and D. L. Turcotte, “Forest fires: an example of self-organized critical behavior,” Science, vol. 281, no. 5384, pp. 1840–1842, 1998. [32] D. L. Turcotte, B. D. Malamud, F. Guzzetti, and P. Reichenbach, “Self-organization, the cascade model, and natural hazards,” Proceedings of the National Academy of Sciences, vol. 99, no. suppl 1, pp. 2530–2537, 2002. [33] L. Palmieri and H. J. Jensen, “The forest fire model: the subtleties of criticality and scale invariance,” Frontiers in Physics, vol. 8, p. 257, 2020. [34] J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models. Cambridge University Press, 1999. [35] M. Henkel, H. Hinrichsen, and S. Lübeck, Non-Equilibrium Phase Transitions. Volume 1: Absorbing Phase Transitions. Springer, 2009. [36] A. Clauset, C. R. Shalizi, and M. E. Newman, “Power-law distributions in empirical data,” SIAM review, vol. 51, no. 4, pp. 661–703, 2009. [37] M. Marsili, I. Mastromatteo, and Y. Roudi, “On sampling and modeling complex systems,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2013, no. 09, p. P09003, 2013. [38] M. Gerlach and E. G. Altmann, “Testing statistical laws in complex systems,” Physical review letters, vol. 122, no. 16, p. 168301, 2019. [39] M. Serafino, G. Cimini, A. Maritan, A. Rinaldo, S. Suweis, J. R. Banavar, and G. Caldarelli, “True scale-free networks hidden by finite size effects,” Proceedings of the National Academy of Sciences, vol. 118, no. 2, p. e2013825118, 2021. [40] A. Pikovsky, M. Rosenblum, and J. Kurths, “Synchronization: a universal concept in nonlinear science,” 2002. [41] K. G. Wilson, “The renormalization group and critical phenomena,” Reviews of Modern Physics, vol. 55, no. 3, p. 583, 1983. [42] V. Loreto, L. Pietronero, A. Vespignani, and S. Zapperi, “Renormalization group approach to the critical behavior of the forest-fire model,” Physical review letters, vol. 75, no. 3, p. 465, 1995. [43] N. Goldenfeld, Lectures on phase transitions and the renormalization group. CRC Press, 2018. [44] G. Jona-Lasinio, “The renormalization group: A probabilistic view,” Il Nuovo Cimento B (1971-1996), vol. 26, no. 1, pp. 99–119, 1975. [45] P. Bak, K. Chen, and C. Tang, “A forest-fire model and some thoughts on turbulence,” Physics letters A, vol. 147, no. 5-6, pp. 297–300, 1990. [46] P. Bak, How nature works: the science of self-organized criticality. Copernicus Books, New York, 1996. [47] D. Rybski, V. Butsic, and J. W. Kantelhardt, “Self-organized multistability in the forest fire model,” Phys. Rev. E, vol. 104, p. L012201, Jul 2021. [48] S. Pueyo, “Self-organised criticality and the response of wildland fires to climate change,” Climatic Change, vol. 82, no. 1, pp. 131–161, 2007. [49] A. Staal, E. H. van Nes, S. Hantson, M. Holmgren, S. C. Dekker, S. Pueyo, C. Xu, and M. Scheffer, “Resilience of tropical tree cover: The roles of climate, fire, and herbivory,” Global Change Biology, vol. 24, no. 11, pp. 5096–5109, 2018. [50] R. D. Zinck, M. Pascual, and V. Grimm, “Understanding shifts in wildfire regimes as emergent threshold phenomena.,” The American Naturalist, vol. 178, no. 6, pp. E149–E161, 2011. PMID: 22089877. [51] M. De Domenico, V. Nicosia, A. Arenas, and V. Latora, “Structural reducibility of multilayer networks,” Nature communications, vol. 6, no. 1, pp. 1–9, 2015. [52] M. De Domenico, C. Granell, M. A. Porter, and A. Arenas, “The physics of spreading processes in multilayer networks,” Nature Physics, vol. 12, no. 10, pp. 901–906, 2016. [53] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, “Multilayer networks,” Journal of complex networks, vol. 2, no. 3, pp. 203–271, 2014. [54] T. E. Harris, “Contact interactions on a lattice,” Annals of Probability, vol. 2, no. 6, 1974. [55] J. Marro and R. Dickman, Nonequilibrium phase transitions in lattice models. Cambridge University Press, 2005. [56] D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” The journal of physical chemistry, vol. 81, no. 25, pp. 2340–2361, 1977. [57] P. V. Martı́n, V. Domı́nguez-Garcı́a, and M. A. Muñoz, “Intermittent percolation and the scale-free distribution of vegetation clusters,” New Journal of Physics, vol. 22, no. 8, p. 083014, 2020. [58] J. E. Keeley, J. G. Pausas, P. W. Rundel, W. J. Bond, and R. A. Bradstock, “Fire as an evolutionary pressure shaping plant traits,” Trends in Plant Science, vol. 16, pp. 406–411, Aug. 2011. [59] P. Jyoteeshkumar reddy, S. E. Perkins-Kirkpatrick, and J. J. Sharples, “Intensifying australian heatwave trends and their sensitivity to observational data,” Earth’s Future, vol. 9, no. 4, p. e2020EF001924, 2021. [60] T. A. Fairman, C. R. Nitschke, L. T. Bennett, T. A. Fairman, C. R. Nitschke, and L. T. Bennett, “Too much, too soon? A review of the effects of increasing wildfire frequency on tree mortality and regeneration in temperate eucalypt forests,” International Journal of Wildland Fire, vol. 25, pp. 831–848, Sept. 2015. [61] Z. L. Steel, D. Foster, M. Coppoletta, J. M. Lydersen, S. L. Stephens, A. Paudel, S. H. Markwith, K. Merriam, and B. M. Collins, “Ecological resilience and vegetation transition in the face of two successive large wildfires,” Journal of Ecology, 2021. [62] R. H. Nolan, L. Collins, A. Leigh, M. K. Ooi, T. J. Curran, T. A. Fairman, V. R. de Dios, and R. Bradstock, “Limits to post-fire vegetation recovery under climate change,” Plant, cell & environment, 2021. [63] V. de L. Dantas, M. A. Batalha, and J. G. Pausas, “Fire drives functional thresholds on the savanna–forest transition,” Ecology, vol. 94, no. 11, pp. 2454–2463, 2013. 5 [64] I. Oliveras and Y. Malhi, “Many shades of green: the dynamic tropical forest–savannah transition zones,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 371, no. 1703, p. 20150308, 2016. [65] K. Cranmer, J. Brehmer, and G. Louppe, “The frontier of simulation-based inference,” Proceedings of the National Academy of Sciences, vol. 117, no. 48, pp. 30055–30062, 2020. [66] G. Nicoletti and D. M. Busiello, “Mutual information disentangles interactions from changing environments,” Physical Review Letters, vol. 127, no. 22, p. 228301, 2021. [67] G. Nicoletti, S. Suweis, and A. Maritan, “Scaling and criticality in a phenomenological renormalization group,” Phys. Rev. Research, vol. 2, p. 023144, May 2020. [68] G. Nicoletti, “Stochastic simulation of a spatial forest-fire model,” Zenodo, Jan. 2023. [69] L. Saravia, “Code for download and Data for Burned Area Australia,” Zenodo, Jan. 2023. [70] R Core Team, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2021. [71] S. M. Bhattacharjee and F. Seno, “A measure of data collapse for scaling,” Journal of Physics A: Mathematical and General, vol. 34, no. 33, p. 6375, 2001.