Intermittent percolation and the scale-free distribution of vegetation clusters

Understanding the causes and effects of spatial patterns of vegetation is a fundamental problem in ecology, especially because these can be used as early predictors of catastrophic shifts such as desertification processes. Empirical studies of the vegetation cover of drylands and semiarid regions have revealed the existence of vegetation patches of broadly diverse sizes, such that the probability distribution of their sizes can be fitted by a power law, i.e. vegetation patches are approximately scale free. Different explanatory mechanisms, such as plant-plant interactions and plant-water feedback loops have been proposed to justify such a finding, yet a full understanding has not been reached. Using a simple individual-based model for vegetation dynamics, we show that environmental variability --a well-recognized characteristic feature of semiarid environments-- promotes the emergence of power-law vegetation patches in a robust way (i.e. for a broad range of parameter values). Furthermore, these are related to a percolation phenomenon that occurs in an intermittent or fluctuating way. The model also reveals that the emerging power-law exponents depend on the overall vegetation-cover density exhibiting the same trend as empirical data and, remarkably, the fitted exponents agree reasonably well with the ones measured in well-known empirical studies. We conclude that environmental variability plays a key role in the formation of such vegetation patterns. From a practical viewpoint, this may be of importance to predict the effects that changes in environmental conditions may have in real ecosystems. From a theoretical side our study shed new light on intermittent percolation phenomena occurring under fluctuating conditions.


Introduction
Living organisms are typically distributed in space neither uniformly nor randomly, but forming structured patterns of aggregation. Understanding the causes and effects of such spatial aggregation is one of the most fundamental problems in ecology [25,12,23,49,50]. An archetypical example is found in arid and semiarid environments, which are characterized by the presence of vegetation patches forming either regular or irregular patterns even in seemingly spatially homogeneous environments. [42,22,47,18,24]. A large number of hypotheses have been proposed to explain the origin of such patterns. These include, among others, the interplay between competition and facilitation mechanisms in plant-plant interactions and plant-water feedbacks [45,19,26,27,28], but, still, we are far from reaching arXiv:1912.08369v1 [cond-mat.stat-mech] 18 Dec 2019 a full understanding. Importantly, variations in the shape of vegetation patterns as well as in the density cover have been proposed as early-warning indicators of catastrophic shifts, such as desertification [46,43,19,20,4]. Thus, understanding, characterizing, and eventually controlling vegetation patterns is a problem of outmost theoretical and practical relevance.
While some of the patterns formed by vegetation, such as stripes, have a characteristic size scale, ground-breaking empirical studies showed a decade ago that vegetation aggregates in semi-arid environments lack a well-defined characteristic scale (i.e. vegetation patches of greatly diverse sizes can be found together) [45,19]. More specifically, in this latter situation patch-size distributions can be fitted by power-laws, the hallmark of scale-free systems [34,30,51,32]. From a theoretical perspective, scale-free cluster-size distributions are often the fingerprint of critical points [5,16,54,9,15]. A classical and archetypical example of a critical point appears in the problem of standard, also called "isotropic", percolation.
In the simplest percolation setup there are two possible states for the sites of a given lattice (which represents the physical space): empty or occupied. In the percolating phase the largest aggregate of occupied sites spans the whole system; on the other hand, in the non-percolating phase only localized finite clusters exist. In between these two phases -right at the percolation threshold or critical point-clusters of all possible sizes -from tiny to system-size-spanning ones-emerge, generating a power-law (scale-free) distribution of cluster sizes with well-known universal exponent values and characteristic scaling features [54,9].
If the ecological processes giving rise to patterned vegetation patches in semiarid ecosystems would follow dynamical rules operating close to a percolation phase transition [44,36], scale invariance should be universal; in other words, there should emerge robust (detail-independent) power-law exponents [5,16]. However, contrary to this expectation, the empirically found distributions of patch sizes can be fitted by power laws whose exponents depend on the level of vegetation coverage in the considered area [45,19,4]. This variability in scaling exponents suggests a lack of universality and sheds doubts on the -maybe exceedingly naive-alleged connection with standard/isotropic percolation and universality, calling for more elaborated theoretical approaches to understand the dynamics and patterns of these ecological communities.
Diverse relatively-simple theoretical models have been proposed to account for the empirical finding of scale-free vegetation patches [45,19,26,27,28]. Some of these include three possible states to describe the system: occupied, empty, and disturbed/degraded (e.g. unavailable to be colonized by vegetation) sites [19,44]. This is in contrast with the case of standard/isotropic percolation models as described above [54,9]. In such three-state models, an external control parameter (which can be interpreted as the grazing pressure, the yearly rainfall, etc.) regulates the level of vegetation cover under stationary conditions and as it is varied -e.g. to mimic harsher conditions-a phase transition from the vegetated to the deserted state appears [44,19,20]. Remarkably, it was found -relying mostly on computational analyses-that cluster size distributions resembling power-laws appear in some of these models for a wide range of external conditions (control parameter values) rather than just around a fixed percolation or critical threshold [44,36]. Thereby, the influential concept of "robust criticality" was introduced by Pascual and colleagues to describe this phenomenon [44,36]. Some other models with only two states (occupied and empty) and some type of facilitation mechanism -as, for instance, a positive feedback between vegetation and water supplies-were claimed to support the idea of robust criticality [45,27]. From a theoretical side it is still not clear where such a generic or robust scaling comes from, nor if it is just a transient (not asymptotic) effect. In fact, van den Berg et al. shed some doubts on the theoretical foundations of the concept by demonstrating analytically that in some simple version of these models scaling occurs at a single point, i.e. there is an underlying sharp transition and not a broad range of criticality [55,56].
The aim of the present work is therefore to scrutinize whether some type of generic scaling -occurring in broad regions of parameter space-can actually emerge in very simple dynamical models (such as the contact process [16]) as well as to uncover the basic mechanisms promoting the emergence of scale invariance in the vegetation-patch sizes within broad parameter ranges.
To this end, we introduce a simple spatially-explicit individual-based model for vegetation dynamics and analyze its behaviour. In particular, we examine if the emergence of scale-free vegetation patches requires a fine tuning of the control parameter (i.e. if it behaves as a standard percolation-like transition) or if, on the contrary, scale-free size distributions emerge in a robust way for a broad range of parameter values. We show that implementing temporal environmental variability in the model -an ingredient inherent to actual ecological systems and of special relevance in semi-arid environments-induces the emergence of scale-free patch-size distributions in a wide region in parameter space , while such a region shrinks to a single point in the absence of temporal variability. Thus, generic scaling in this case emerges as a sort of intermittent or fluctuating percolation phenomenon, in which clusters nearby a percolation threshold fluctuate in size owing to external variability. Remarkably, we show that our model -even if with some limitations induced by its simplicity-is able to reproduce with reasonable accuracy the values and overall dependence on vegetation density of some empirically observed non-universal scaling exponents.

Model building
The parsimonious individual-based model that we analyze is defined on a two-dimensional square lattice, with each lattice site being in one of two possible states: empty (0) or occupied (1). At each discrete time step an occupied site is randomly selected and, with probability p, a new offspring is created at a random nearest-neighbor provided it was empty [29,16]. Alternatively, with probability 1 − p the selected occupied site is emptied; i.e. an individual dies. The time counter is increased as t → t + ∆t with ∆t = 1/N (t), where N (t) is the total number of occupied sites at time t. This way, we ensure that in a (Monte Carlo) step all are updated once on average. This stochastic rules are nothing but the "contact process", a prototypical model that captures the essence of a large variety of spreading systems in physics, biology, epidemiology, [29,16,15] as well as in spatial ecology [12,33,13,61,60,7]. It can be summarized as the following set of processes: where A indicates "particles" or "plants" in our interpretation.
In the special case of a well-mixed population in which all sites are taken to be neighbors of all others and assuming an infinitely large system (i.e. in the mean-field case) the dynamics of the overall density of occupied sites ρ(t) = 1/N N i s i obeys [5,16,15]: in some unspecified time units 1 . This equation exhibits a transcritical bifurcation at the reproduction probability p c = 1/2. This is a directed-percolation critical point [29,16,35,15], above which the population survives for arbitrarily long time in sufficiently large lattices reaching a stationary density ρ s = (2p − 1)/p. Below such a value, the only stable solution is the empty state ρ s = 0. Thus the critical point separates an empty (or "absorbing") phase from an occupied (or "active") phase.
This type of exact solution cannot be straightforwardly extended to analytically solve the case of a two-dimensional lattice in simple terms; however, it is well-known that the overall phenomenology is qualitatively (though not quantitatively) preserved [29,16]. This is illustrated in Figure 1 which shows the phase diagram for the two-dimensional contact process resulting from computer simulations. Observe, in particular, that the directed-percolation critical point -separating the empty from the active phase-is shifted to p c ≈ 0.6225 owing to fluctuations, i.e. to the finite size of local neighborhoods [29,16,15].
Most of the literature on the contact process is focused in analyzing the non-trivial features of this absorbing-active phase transition [29,16,35,15]. However, here we are interested in another "hidden" critical point whose very existence is usually neglected. This -as illustrated schematically in Figure 1-is a percolation-like transition at which a cluster of occupied sites percolates in space [55,56]. More specifically, this second transition is -as we shall showa standard/isotropic percolation phase transition [54,9] above which, at any given time within the steady state, the probability to find a percolating cluster of adjacent occupied sites spanning the whole two-dimensional space is larger than zero. In two-dimensions, one needs to reach stationary activity densities of the order of ρ ≈ 0.56 (deep inside the active phase, i.e. p ≈ 0.725 > p c ) to observe the emergence of a spatially percolating cluster ( Fig. 1). Let us stress, even at the risk of being redundant, that this transition at which clusters percolate in space is conceptually different from the absorbing-to-active directed-percolation transition, at which clusters of occupied sites can percolate in time.
We first study the nature of the percolation transition within the active phase of the "pure" contact process. Although our main focus will be on a variant of the contact process aimed at modelling environmental variability, which we proceed now to describe.
In the temporally-disordered contact process, at each Monte Carlo timestep the reproduction probability, p(t), takes a uniformy-distributed random value in the interval [p + σ,p − σ]. The values ofp and σ are kept constant in time in each simulation, butp can change for different realizations to explore diverse regions of parameter space, and in particular to sample different densities.p is taken within the interval [σ, 1 − σ] so that p(t) ∈ [0, 1]. The width σ is kept as a control parameter characterizing the environment variability, e.g. the year-to-year demographic (i.e. birth and death) probabilities 2 . In this temporally-disordered version of the model (fixing σ = 0.25), the absorbing-to-active σ=0.00 σ=0.25 (separating empty from populated states), the stationary density is ρ = 0, and it is not possible to measure any non-trivial steady-state spatial cluster distribution. Observe that within the active phase -in which activity propagates, reverberating indefinitely across time-clusters of adjacent active sites may (i) not percolate in space (as occurs, e.g. for p = 0.65, which corresponds to ρ 0.26; see lower inset) or (ii) do so, deeper into the active phase (as occurs for e.g. p = 0.725, corresponding to ρ 0.56 for large system sizes; see upper inset). In the insets we have highlighted the largest cluster (dark orange color) in snapshots for the above two cases, illustrating the absence of percolation for ρ 0.26, and the existence of a percolating cluster for ρ 0.56. For the temporally heterogeneous case, the displayed density is the average over time for an average reproduction probabilityp ∈ [σ, 1 − σ], such that p(t) =p + σ lay in [0, 1]. In this case, the absorbing-to-active critical point is shifted to a value (p c ≈ 0.6350) larger than that of the pure model. critical point is located at p c ≈ 0.6350, a value slightly larger than the corresponding value for the homogeneous (or "pure") case p c ≈ 0.6225 (see Fig.1) [59].
Let us stress that in order to properly compare our computational results for either of the considered models with field observations for vegetation cover at a given vegetation density ρ [45,19], we need to consider the "constant − ρ ensemble". This means that statistics should be collected not for a constant birth probability p (which is difficult to measure empirically) but for constant density ρ (a more accesible quantity). Although in the infinitely-large system-size limit, "constant-p ensemble" and "constant − ρ ensemble" are expected to be equivalent [17,16], in finite systems -as there is not a one-to-one map between p and ρresults can be (and actually are) different in both ensembles.
In practical terms, to obtain results for a given density ρ, we run computer simulations for different values of p and for long-enough times as to reach the stationary state (typically t > 10 3 ). Unless specifically mentioned, we consider a system size N = 1024 × 1024 which is enough to collect good statistics at a reasonable computational cost. Then, we take a "snapshot" of the final configuration and collect the statistics of patches sizes conditioned to the value of its overall density ρ (using a binning of size ∆ρ = 0.001). As the output, we determine the path-size distribution for each possible density value.

Results
The (isotropic) percolation transition in the "pure" contact process First, we analyze the percolation transition of the homogeneous or "pure" contact process -without additional heterogeneity-and use it as a benchmark for further analyses. In two-dimensional percolation on a square lattice the percolation threshold is p occ c ≈ 0.5927... and at such a critical point the cluster-size distribution decays as with an exponent value τ = 187/91 ≈ 2.05 [54,9]. Importantly, the value of p c can change if other types of latices (e.g. hexagonal) are considered, but the exponent remains unchanged, i.e. it is universal.
Our dynamical models are (obviously) not identical to the standard/isotropic percolation model. Among other aspects, the models here are dynamical and the states of nearby sites are clearly correlated. The question we address now is whether such differences are relevant, i.e. whether the universal aspects of the percolation transition are preserved or not.
To provide an answer to this question we resort to computational analyses. Fig.2A and B show the cluster-size distribution P (s) for a wide range of occupation-densities ρ of the pure contact process (low (A) and high density (B) curves are separated for the sake of visual clarity). Observe that at a critical occupation-density value ρ c ≈ 0.56 (and only at this value) the cluster-size distribution decays as a power-law, i.e. cluster vegetation patches of broadly different sizes emerge. Furthermore, the scaling behavior is fully compatible with the expected universal behavior of standard/isotropic percolation, fitted value for the cluster-size distribution exponent τ is in agreement with the above-mentioned value τ ≈ 2.05 [54,9]. is characterized by a cluster-size distribution P (s) such that P (s) ∝ s −τ with τ = 187/91 ≈ 2.05 and the scaling relationship χ = |p occ − p occ c | −γF (|p occ − p occ c |L 1/ν with ν = 4/3 and γ = 43/18 [54,9]. We show P (s) for low (A) and high (B) densities (plotted separately for the sake of clarity) as obtained in computer simulations of the contact process on a square lattice, for different occupation densities ρ in the range ρ ∈ [0.02, 0.70] (color coded). Panel C shows the curve collapse for four different system sizes, using the scaling exponents of standard/isotropic percolation. We conclude that the contact process displays clean scaling behavior for a critical density ρ c ≈ 0.56, fully compatible with that of the standard/isotropic percolation class.
Other quantities usually employed in studies of standard/isotropic percolation can also be determined at the percolation transition of the pure contact process. In particular, right at the percolation threshold, the mean cluster size χ, defined as where i runs over all clusters, s i is the corresponding cluster size and n occ is the total number of occupied sites; is governed by the general scaling law where ∆ = |p occ − p occ c | is the distance to the critical occupation probability p occ c , L the system linear size, and ν = 4/3 and γ = 43/18 are exactly-known critical exponents [54,9]. This implies that there exists a master curveF(∆L 1/ν ) in which data should collapse for different system sizes L and ∆ values (if the right exponent values are used for ν and γ). Figure 2C shows such a curve collapse for the contact process at fixed density values -at the critical point ρ c ≈ 0.56providing us with strong evidence that it exhibits the same universal scaling behavior as standard/isotropic percolation.
In order to rationalize this conclusion from a theoretical viewpoint one can argue as follows. At the percolation critical point of the contact process -which lies well inside the active phase-occupied sites are not randomly placed, as dynamically-generated spatial correlations exist. However, as the dynamics at such a density is far away from the absorbing-active phase transition (where long-range correlations emerge), correlations are only short-ranged. As a consequence and from a renormalization-group perspective, at a sufficiently "coarse-grained" scale, the effective density becomes uncorrelated. Thus, standard/isotropic percolation scaling emerges at sufficiently large scales. In other words, the presence of weak short-range correlations changes the location of the critical density but does not affect the large-scale behavior of percolation, in perfect agreement with what computationally observed.
Summing up, our simple pure model exhibits scale invariance in its cluster-size distribution if, and only if, it is fine tuned to have a density close to the percolation threshold. Around such a density it exhibits scaling features of standard percolation, and no trace of generic scale invariance for arbitrary density values, for which exponential behavior (with a characteristic scale) is observed.
Constant-ρ ensemble properties with temporal heterogeneity.
When temporal heterogeneity, σ > 0, is introduced, the scenario described in the previous section for the case of the pure contact process (σ = 0) changes dramatically, as we show in what follows.
Unless explicitly said, all simulations are ran for σ = 0.25, without loss of generality. Given the subtleties usually encountered when determining whether a given distribution follows a power law or not, let us stress that in the forthcoming analyses we employ the standard most-stringent statistical tests -following the methods introduced in [10] (see also Appendix)-to fit power laws and to decide whether they provide better fits than e.g. exponentials.
In particular, for each distribution P (s) at a given density ρ we obtain the parameters C, τ and ∆s in Eq.(3) that maximize the normalized log-likelihood of the fit (see Appendix). Our analyses reveal that cluster-size distributions conform to a power-law curve for at least two (and up to four) decades (see Fig. 3A and B, and Table 1). In this case, contrary to the pure contact process, the scale-free behavior of the cluster-sizes holds within a wide range of densities, with continuously varying exponents. Simple eye inspection of Fig. 3A and B already reveals that the exponent τ , characterizing the observed power-law decay, displays a non-monotonic variation with population density ρ (see Fig.  3C). In particular, the fitted value decreases from τ 2.59 (for ρ = 0.02) up to τ 1.63 (for ρ ≈ 0.5), that value τ increases again up to τ 3.12 (for ρ ≈ 0.7) (see Fig. 3C and Table 1).
Observe that -similarly to the standard/isotropic percolation case-the average cluster size of finite clusters (i.e. the mean of these distributions) increases monotonously as the control parameter is increased in the subcritical phase and it starts decreasing once the percolation threshold is reached and finite clusters start merging together or with the percolating giant cluster (if it exists) [9]. However, it is remarkable that opposed to the standard/isotropic percolation case, power laws provide a good fit for the cluster-size distribution not only at the critical point but also in a wide parameter range above and below such a value, i.e. we observe scale-free cluster sizes through a wide range of densities.
Owing to this peculiarity, the percolation threshold ρ c cannot be straightforwardly obtained -as usually done-just by looking for scale-free cluster-size distributions, since there is not a unique density where scaling is observed. Thus, we resort to obtain an estimate of the critical point given by the density for which the fit to a power law extends for the largest possible range of cluster sizes with a fit exhibiting a maximal normalized log-likelihood. This is equivalent to looking for the point with largest averaged finite-cluster size, as usually done in standard/isotropic percolation. Such an optimal value turns out to be around ρ ≈ 0.56, roughly speaking, the same occupation density obtained for the isotropic percolation transition of the pure contact process without temporal disorder (see Table I where details of the fits are given).
The exponent fitted at ρ = 0.56 is τ 1.81, slightly smaller than the expected value in two-dimensional percolation τ ≈ 2.05 (see Table 1). Moreover, considering ρ = 0.56 as the percolation critical point, we obtain a scaling collapse (see Fig. 3D) similar to that of the homogeneous case, represented in Fig.2C but with exponents ν 3.5 and γ 4 that are considerably larger than the percolation values (ν = 4/3 and γ = 43/18), suggesting a different type of critical behavior (or universality class). A more precise determination of the critical point and the associated critical exponents is beyond the scope of the present manuscript and will be reported elsewhere. Observe that contrarily to the pure contact process case, power law behaviour can be fitted for a wider range of densities (at least two or more decades for almost all densities, and up to four decades around the critical point). For each density ρ, we obtain a power-law exponent τ such that it maximizes the normalized log-likelihood (see Appendix, [10], panel C) giving raise to a non-monotonic dependence. The best fit is obtained for ρ = 0.56 with an exponent τ 1.81, slightly smaller than the expected for the standard percolation τ = 187/91 2.05. Taking this value as the critical percolation point ρ occ c = 0.56 we obtain a rather good scaling collapse (panel D) for χ = ∆ −γχ (∆L 1/ν ) for values ν 3.5 and γ 4, much larger than expected for standard/isotropic percolation ν = 4/3 and γ = 43/18.

Individual snapshot properties with temporal heterogeneity.
We have just shown that our model equipped with temporal disorder is able to reproduce scale-free cluster-size distributions without the need of careful fine tuning to a threshold point, i.e. clusters of largely different scales emerge in a rather robust way for a broad range of parameter values. Nevertheless, in order to compare with empirical observations, an important caveat remains. The results obtained so far correspond to ensemble averages where many different realizations of the system are considered and averages are performed over many different times in the steady state, with the only restriction that the value of the control parameterp is kept constant. On the other hand, empirically measured exponents of vegetation-patch-size distributions correspond to single snapshots of the (unknown) ecological dynamics.
To bridge this gap, we also computed the statistics of individual (or "instantaneous") snapshots within our model. In other words, we performed a set of simulations of the model for diverse values ofp and took instantaneous pictures of the system configuration in the steady state. For each of such snapshots we compute its overall density and measure the cluster-size distribution (histogram). Obviously, these distributions -not being ensemble averaged-are much noisier than those reported above. Eye inspection of Figure 4 already reveals that the cluster-size distribution of an individual snapshot -even for a density far away from the percolation threshold, ρ = 0.3 in this case-is broadly distributed and can be fitted as a power-law (straight line in the double-logarithmic plot). This panel shows the cluster-size distribution as obtained from panel (A); the distribution is obviously noisy given the relatively small system size, but it is clearly a broad one, and can be fitted as a power law (black dashed line).
A collection of individual-snapshot cluster-size distributions together with the averaged distribution are shown in Fig.5 for different densities. Observe that in some cases the averaged distribution extends to sizes larger than those of individual snapshots. This is just because, for the sake of clarity, not all individual snapshots in the average are plotted. This figure illustrates that even for a fixed density ρ there is a large variability in the histograms of cluster-sizes and that single snapshots may differ notably from the ensemble-averaged (constant-ρ) distribution. In particular, fast exponential-like and slow power-law-like decays can be identified for individual snapshots. Although exponential decays are seen, especially for low densities, power-laws are present for snapshots of all density values; however, they are more prominent and broad in the neighborhood of the percolation transition (see Table II where details of the fits are reported). In this case of histograms for individual snapshots, we also compared the goodness of power-law fits with that of exponential-distribution fits P (s) = Ce −αs . For most of the snapshots power-law distributions fit better the data than exponential ones, i.e. power-laws give a higher normalized log-likelihood in the majority of the cases, but not always, especially for low densities ρ ∈ [0.04, 0.24] (see fourth column of Table 2). However, even for densities in this interval, power-law fits are better in at least 60% of the snapshots in all cases (see Table 2). On the other hand, for intermediate densities, ρ ∈ [0.26, 0.68], power laws outperform exponential fits in almost all snapshots. Observe that in spite of the limitations imposed, the finite system size and the lack of an averaging procedure, power-laws are observed for at least 1.5 decades in most of the cases (fifth column of Table 2).
The average of the exponents obtained by fitting power-law distributions for individual snapshots (obtained excluding samples that are best fitted by exponentials) exhibits a non-monotonic behavior as a function of the vegetation density (see gray line of Fig. 6). The mean exponent τ (ρ) reaches a minimum value τ 1.54 < τ c at a density value below the percolation threshold ρ 0.44 < ρ c , and increases up to τ 4 for larger and smaller densities.
We also performed simulations for smaller values of the noise amplitude σ, for which the shape of the curves remain similar and the widths of the distributions of τ -values (i.e. the errorbars in Fig.6) decrease. Actually, in the limit in which σ goes to 0 one should recover power-law scaling only in a vanishing region around the critical point. However, a detailed comparison between the resulting curves for diverse values of σ is not straightforward given that -in the way we chose to parametrize-the possible values ofp are limited by the choice of σ and thus the sampling of values of ρ is also affected by σ. Disentangling these effects from the final results is therefore not straightforward in the present setting; thus, we do not pursue further such a goal here.
Comparison with empirical data.
We now compare the obtained exponents τ of the snapshots distributions P (s) with available empirical data [45,19]. In particular, the black dots of Fig. 6 correspond to empirical data for different vegetation-cover densities. Remarkably, field data are clearly compatible with our model results (for large enough environmental variability σ = 0.25) and within one standard deviation (p-value p = 0.04 < 0.05 of Pearson's χ 2 statistics).
In the light of these results, we can conclude that a very simple model of vegetation dynamics -when equipped with temporal environmental variability as a key ingredient-is able to produce highly structured (vegetation) configurations which can be well fitted by power laws in a rather robust though non-universal way. Furthermore, and rather remarkably, the model also shows exponents that change with the cover density, having values compatible with those found in empirical communities.

Discussion
The study of vegetation patterns have a long tradition in ecology and environmental sciences. This interest has been reinforced in recent years owing to the growing awareness that the statistics of such patterns can be used as a predictor of early warning of ecological transitions or regime shifts such as desertification in semiarid environments. From the empirical side, two very influential papers reported a decade ago that vegetation patterns in some arid and semiarid ecosystems are scale free, thus lacking characteristic scales [45,19]. Puzzled by these field observations, theoreticians were challenged to understand the origin of such scale-free patterns and alternative explanatory mechanisms have been proposed to account for them. For instance, some models put the emphasis on (i) existing plant-plan interactions (competition and facilitation), (ii) feedback loops between vegetation dynamics and abiotic resources such as water [45,26], or (iii) the presence of heterogeneity in the environment [19,48,27,28].
For the sake of generality, let us stress here that the theoretical challenge of modelling approaches is to explain how scale-free patterns can possibly emerge in a robust way, i.e. without the need of invoking a high degree of parameter  Figure 6: Fitted exponents agree with field observations. We fitted power-law regimes of cluster-size distributions P (s) of snapshots for different densities ρ and σ = 0.25 (see Appendix and Table 2 for details). We obtained the mean (gray lines) and standard deviation (shaded area) of the exponents τ obtained for up to 200 snapshots (see sixth columns of Table 2) for density values [0.04, 0.7] equally spaced in increments of size 0.01. Observe that empirical data [45,19] (black dots) are in good agreement with the model predictions (at least, for enough environmental variability σ = 0.25).
fine tuning. In particular, experience in statistical physics tells us that scale-free patterns typically emerge at very special points of the parameter space, such as critical points. To circumvent this theoretical conumdrum and in order to explain generic power-law decays, two different types of approaches are considered. The first one consists in considering mechanisms of self-organization, inspired in self-organized criticality [3,40,11,6], which would justify why critical-like scale-free features appear without the apparent need of parameter fine tuning (although they usually require an unrealistic separation of time scales). This type of approach leads, however, to universal -as opposed to continuously varying-exponents. The second strategy consists in devising models in which scale-free patterns -with their concomitant power law distributions-emerge for broad ranges of parameter values [14]. Such type of models already exist in statistical mechanics (for an in-depth discussion of these issues we refer to [14]), but, to the best of our knowledge, they are not well-suited to be translated to the ecological problem under study here.
A particularly relevant idea in this context is that of "robust criticality", stating that scale-free clusters of vegetation are akin to percolation clusters at a percolation phase transition, but that -for reasons to be fully clarified-it could be observed in simple models of vegetation dynamics without the need of fine tuning exactly to the percolation transition [37,38,44,36,39]. In other words, it is claimed that these models exhibit a broad region of power-laws around a transition point. However, the origin of such a broad region of criticality is still controversial in the sense that there is no actual theoretical understanding of why generic scale invariance could possibly emerge.
In this paper we have proposed a very simple model that is able to justify in a clean way the emergence of scale-free vegetation patterns in a broad range of parameter values.
Our model is individual-based and spatially explicit, based on the two-dimensional contact process. This model is the simplest possible to mimic vegetation dynamics (as a birth-and-death process). It has been very deeply studied in the literature and it is well-known to exhibit a quiescent-to-active phase transition, separating a quiescent phase on which all activity/vegetation ceases/disappear, from another in which activity/vegetation survives indefinitely with a given averaged density. In other words, the model exhibits a "directed-percolation" type of phase transition separating phases where activity may, or may not, percolate in time.
On the other hand, it is often neglected that there is another phase transition hidden in this prototypical model. It occurs within the active phase (as illustrated in Figure 1) and separates a phase in which the largest cluster of active sites percolates in space from another one in which it does not. We have confirmed computationally that, indeed, as theoretically predicted by van den Berg et al. [55,56] the transition is sharp, i.e. there is not such a thing as a broad range of scaling for such a simple model. Moreover, we have shown for the first time that such a transition is perfectly described by the standard scaling theory of (isotropic) percolation, with its well-known set of scaling exponents. Let us emphasize that such a result is not trivial as in standard percolation sites are occupied in a completely random way, while in the active phase of the contact process activity occupied sites are spatially correlated. Readers familiar with the renormalization-group theory could have correctly anticipated though -employing renormalization-group arguments-that such short-ranged correlations should not be able to affect the scaling exponents. Thus, summing up, the standard contact process only exhibit scale-free clusters at a very precise point in parameter space and is therefore not an adequate model to explain the broad/generic emergence of scale-free vegetation patterns in a robust way.
Inspired by previous work of our group and others showing that temporal variability can lead to generic power laws in some simple models [59,64], we moved on to study a variant of the contact process implementing temporal variability in external conditions, i.e. a randomly changing control parameter. With this simple -and ecologically well justifiednovel ingredient the model changes dramatically.
i) First of all, the percolation phase transition within the active phase is not in the standard/isotropic percolation universality class, i.e. its scaling exponents are significantly different form the standard ones (see below for a further discussion of this fact).
ii) Second, generic scale-free distributed clusters emerge for a broad range of parameter values and not just at the critical point.
iii) Third, individual snaphots of the system state -i.e. single configurations in which no averaging is performedalso exhibit clusters of highly variable sizes that can be fit as power laws for a broad range of possible vegetation densities. Even more remarkably, the obtained power-law exponents are density dependent -as also occur in empirical observations-and are in reasonably good agreement with those observed in the field for each single density for which empirical results were accessible to us.
Let us now briefly discuss each of these three results.
i) To understand in a qualitative way why the scaling at the percolation point does not coincide with that of standard/isotropic percolation, one could argue as follows. When unfavorable external conditions appear, small clusters are very likely to be removed, while large ones, although also affected by poor conditions, can remain connected or divided in smaller-sized clusters. When favorable periods return, the latter clusters become more compact and/or are re-merged with others. Contrarily, the removed clusters are not able to re-appear from the empty sites, so that empty regions persist. As a result of these effects, the active sites in the active phase are strongly correlated and not equivalent to a standard/isotropic percolation process. Of course, in order to turn these hand-waving arguments into a more formal proof, one would need to perform some type of renormalization-group calculation, which is far from our scope here. However, let us remark, that usually the presence of temporal disorder dramatically affects the universal behavior of well known universality classes [63,1] (including the quiescent-active critical point of the contact process, which is dramatically affected by temporal noise [31,8,59]). Thus, it does not come as a surprise that our model equipped with temporal disorder differs essentially from standard percolation. In any case, a systematic account of this (novel?) critical universality class, establishing all critical exponents and scaling relationships with good accuracy and precision, is left for a future publication.
ii) The emergence of a broad region of densities for which scale-free behaviour emerges is in line with the idea of generic scaling emerging in the presence of temporal variability. Actually, it was recently shown that in the presence of temporal heterogeneity a novel type of phase may emerge in systems with an absorbing-active phase transition such as the contact process. Such a novel phase was named temporal Griffiths phase [59] in analogy with the well-known Griffiths phases that emerge in the presence of spatial (quenched) heterogeneity [62]. In particular, it was shown in [59] (see also [58,64]) that for models similar to the contact process studied here, different quantities scale as power laws, not just at a critical transition point, but in a whole broad region in parameter space in the presence of temporal disorder.
iii) Directly related to our findings is the important work of Manor and Shnerb who discussed -in a mean-field set up different from our spatially explicit one-how multiplicative (environmental) noise can lead to scale-free distributed cluster sizes [26,28]. In a nutshell, each existing cluster experiences a multiplicative random-walk process in which the change in cluster size either increases or decreases at each time step, but, crucially, it does so in an amount which is proportional to its present size, i.e. in a multiplicative way. It is very well-known that multiplicative processes lead generically to power-law distributions (of sizes in this case) [30,41,53]. We believe that a multiplicative process similar to this one is at the basis of the power-laws found in our model, but it can not be analyzed in simple mean-field terms.
Thus, the reason why our model produces generic power-law distributions is twofold: (i) on the one hand, nearby the standard percolation phase transition clusters of very different sizes are generated and (ii) owing to the multiplicative process induced by changing environmental conditions such cluster fluctuate wildly. The combination of these two features leads to -at least approximate-scale-free cluster distributions. In other words, the model can exhibit a sort of intermittent percolation in which the overall density fluctuates around the critical (isotropic) percolation threshold. This sweeping through the transition (or even near the transition) may suffice to generate power-laws in a rather robust fashion. In fact, D. Sornette proposed years back a generic mechanism like this -alternative to self-organized criticalityto explain the robust emergence of power-laws [52].
Finally, in order to test the robustness of our conclusions, and following some proposals in the literature, we have constructed versions of our temporally-varying model implementing spatial heterogeneity (spatial-dependent p) and a facilitation mechanism (p depends on the number of occupied nearest neighbors). Preliminary studies reveal that results do not change significantly with respect to what is reported here for the pure case. More detailed and careful analyses of such extensions will be presented elsewhere.

Conclusion
The conclusion of the present study is that a very simple individual-based model with fluctuating external conditions is able to explain in a simple but robust way the emergence of scale-free vegetation patterns. The exponents associated with the concomitant power-law distributions are non-universal and density dependent, and agree remarkably well with those fitted to empirical data. This conclusion does not imply that other mechanisms, such as plant-plant interactions or plant-water feedbacks are not important to explain vegetation patterns in general, but reveals that temporal variability is a key ingredient to explain the emergence of scale-free vegetation patterns in semiarid environments.

Acknowledgements
We acknowledge the Spanish Ministry and Agencia Estatal de investigación (AEI) through grant FIS2017-84256-P (European Regional Development Fund), as well as the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), ref. SOMM17/6105/UGR for financial support. We thank P. Moretti, V. Buendía, G. Barrios, and P. Villegas for a critical reading of the manuscript.

Appendix
Fitting power laws We fit cluster-size distributions P (s) -both averaged over different realizations and times and for individual snapshotsas power-laws P (s) = Cs τ within an optimal interval [s min , s max ]. The exponent τ , the normalization constant C as well as the values of s min and s max are fitted by maximizing the normalized log-likelihood assuming Poissonian counts ln L = (1/N ) i [s i ln P (s i ) − P (s i ) + ln(s i !)], where N is the number of non-zero cluster sizes in a given range of s [10,21,2]. To avoid overfitting, we discard snapshots whose optimal range [s min , s max ] includes less than 10 points. In all cases, we compare the power-law fit with an exponential one P (s) = Ce −cs within the obtained interval [s min , s max ]; on the case of snapshots, we report on the number of trials for which the power-law constitutes a better fit than the exponential.  Table 1: Fits for averaged cluster-size distributions P (s) with different densities. Power-laws have been fitted (see Appendix) to averaged P (s) for different densities (first column). We show the fitted exponents (second column), the optimal interval of clusters sizes -s min (third column) to s max (fifth column)-and the corresponding normalized log-likelihood for each density. The best fit is obtained for density ρ = 0.56 while for ρ = 0.48 the power-law fit extends for the broadest range.  Table 2: Fits for individual snapshots with different densities. We fit power-law curves (see Appendix) for snapshots of the model with temporal variability for different densities (as listed in the first column). The second column shows the mean fitted exponents; the third column shows the variance of exponent values (only snapshots whose optimal range considers more than 10 points are considered); the fourth one shows the percentage of snapshots for which power-law fits are better than exponential ones. The fifth column shows the fraction for which the power-law fit extends for at least 1.5 decades.