Intermittent percolation and the scale-free distribution of vegetation clusters

Most species are typically distributed in space neither uniformly nor randomly, but forming patterns of spatial aggregation [1, 2]. Understanding the causes and effects of such spatial aggregation is one of the most fundamental problems in theoretical ecology [3–10]. An archetypical example of vegetation patterns is found in arid and semiarid environments, which are characterized by the presence of vegetation patches, organized in either regular or irregular patterns, even in seemingly spatially homogeneous environments. [11–16]. A large number of mechanisms have been proposed to explain the origin of such vegetation patterns. These include, among others, the interplay between competition and facilitation processes in plant–plant interactions as well as plant-water feedbacks [17–21], but a full understanding of how spatial patterns emerge is still lacking. Importantly, even if this idea remains controversial, 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 [18, 22–25]. Thus, in summary, 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 semiarid environments lack a well-defined characteristic scale (i.e. vegetation patches of greatly diverse sizes can be observed in nearby areas) [17, 18]. More specifically, in this latter situation patch-size distributions can be fitted by power-laws [6, 9], the hallmark of scale-free systems [26–29]. From a theoretical perspective, scale-free cluster-size distributions are often the fingerprint of critical points, at the edge of a phase transition [30–34]. An archetypical example of a critical point, paradigmatic to rationalize scale-free clusters, is 'standard' or 'isotropic' percolation.

In the simplest percolation model, the sites of a given lattice (which represents the physical space) can be in two possible states: either empty or occupied. In the percolating phase, i.e. for large occupation probabilities, the largest aggregate of occupied sites spans the whole system; on the other hand, in the non-percolating phase, i.e. for low occupation probabilities, only finite-size clusters exist. In between these two phases—right at the percolation threshold or critical point—clusters of all possible sizes—from tiny to extremely large ones—emerge, generating a power-law (scale-free) distribution of cluster sizes; this has well-known universal critical exponents and characteristic scaling features [32, 33].

If the ecological processes giving rise to scale-free vegetation patches in semiarid ecosystems would generate—for some unknown reason—clusters such as those close to a percolation phase transition [35, 36], their observed scaling exponents should be universal. In other words, robust power-law exponents, not depending e.g. on vegetation density, should be observed [30, 31]. However, contrary to this expectation, the empirically found distributions of patch sizes can be fitted by power laws whose exponents are system dependent and change with the level of vegetation coverage in the analyzed area [17, 18, 25] (see, however, [2]). This variability in scaling exponents suggests a lack of universality and sheds doubts on the—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 dynamical models have been proposed to account for the empirical finding of scale-free vegetation patches [17–21]. Some of these include three possible states to describe the system at each spatial location: occupied, empty, and disturbed/degraded (i.e. unavailable to be colonized by vegetation) [18, 35]. In such three-state models, an external control parameter (which can be interpreted as the grazing pressure, the yearly rainfall, fire etc.) regulates the level of vegetation cover under stationary conditions. As it is varied—e.g. to model transiently harsher or milder environmental conditions—a phase transition from the vegetated to the deserted state may appear [18, 24, 35]. Remarkably, it has been found—relying mostly on computational analyses—that cluster-size distributions resembling power-laws can be observed in some of these models for a wide range of external conditions (a broad set of parameter values) rather than just around a fixed critical threshold, as in isotropic percolation [35, 36]. Thereby, the influential concept of 'robust criticality' was introduced by Pascual and colleagues to describe this phenomenon [35, 36]. 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 [17, 20]. From a theoretical viewpoint, it is still not clear what is the origin of such generic or robust scaling, nor if it is just a transient (not asymptotic) effect or if, instead, it describes true scale invariance. Actually, van den Berg et al cast doubts on the theoretical foundations of the concept of 'robust criticality' by analytically proving that in a simple version of these models true scaling occurs only at a single point in parameter space, i.e. there is an underlying sharp phase transition and not a broad range of criticality [37, 38]. Thus, understanding the origin of broad (not fine tuned) scaling in this type of systems remains a challenge for theoreticians.

The aim of the present work is therefore to scrutinize whether generic scaling—occurring in broad regions of the parameter space—can actually emerge in very simple dynamical models [31], as well as to help uncovering the basic mechanisms promoting the emergence of scale invariance in vegetation patterns. To this end, we introduce a simple spatially-explicit individual-based model of vegetation dynamics and analyze its behavior. 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 semiarid environments [39]—induces the emergence of broadly distributed cluster sizes (that can be fitted as power laws) in a wide region of the parameter space, while such a region shrinks to a single point in the absence of temporal variability. Hence, 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 approach—even if with some limitations inherent to its simplicity—is able to reproduce the overall dependence on vegetation density of some empirically observed non-universal scaling exponents, thus contributing to shed light on the ecological mechanisms behind such a phenomenology.

1.1. Model building

In two-dimensional percolation on a square lattice the percolation threshold is ${p}_{\text{c}}^{\mathrm{o}\mathrm{c}\mathrm{c}}\approx 0.5927\dots$ and right at such a critical point the cluster-size distribution follows a power law

$$\begin{equation}P\left(s\right)=C{s}^{-\tau }\end{equation} \tag{ 1 }$$

where C is a normalization constant and the exponent takes the exactly known value τ = 187/91 ≈ 2.05; in finite systems this power law is cutoff by a rapidly decaying scaling function [32, 33]. Importantly, the value of ${p}_{\text{c}}^{\mathrm{o}\mathrm{c}\mathrm{c}}$ can change if another type of lattice of local geometry is considered, but the exponent τ remains unchanged, i.e. it is universal.

The individual-based model that we analyze is defined on a two-dimensional square lattice (mimicking the surface of a semiarid landscape) in which each lattice site, j, (i.e. each spatial location) can be in one of two possible states: empty, n(j) = 0, or occupied n(j) = 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 [31, 40]. Alternatively, with probability 1 − p the selected occupied site is emptied (i.e. the 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 sites are updated once, on average. We consider a system size of N = 1024 × 1024, unless otherwise indicated.

These stochastic rules define the 'contact process', a prototypical model that captures the essence of a large variety of spreading systems in physics, biology, epidemiology, [31, 34, 40] as well as in spatial ecology [4, 10, 41–45]. It can be summarized as the following set of processes:

$$\begin{equation}\begin{matrix}\hfill A\hfill & \hfill p{\to }\hfill & \hfill 2A\hfill & \hfill \left(\text{birth}\right),\hfill \\ \hfill A\hfill & \hfill 1-p{\to }\hfill & \hfill 0\hfill & \hfill \left(\text{death}\right)\hfill \end{matrix}\end{equation} \tag{ 2 }$$

where A represents individual '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' limit) the dynamics of the overall density of occupied sites $\rho \left(t\right)={\sum }_{j=1}^{N}n\left(j\right)/N$ obeys

$$\begin{equation}\dot {\rho }\left(t\right)=p\rho \left(t\right)\left(1-\rho \left(t\right)\right)-\left(1-p\right)\rho \left(t\right)=\left(2p-1\right)\rho \left(t\right)-p{\rho }^{2}\left(t\right)\end{equation} \tag{ 3 }$$

in unspecified time units [30, 31, 34].⁵ This mean-field equation exhibits a critical point at the reproduction probability p_c = 1/2 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. Therefore, the critical point separates an empty (or 'absorbing') phase from an occupied (or 'active') phase, i.e. this is a directed-percolation critical point [31, 34, 40, 47].

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 [31, 34, 40, 47]. 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 > 0.5 due to fluctuations, i.e. to the finite size of local neighborhoods.

Zoom In Zoom Out Reset image size

Most of the literature on the contact process is focused on analyzing the non-trivial features of this directed-percolation absorbing-active phase transition [31, 34, 40, 47]. However, here we are interested in another 'hidden' critical point whose very existence is usually ignored. This, as sketched in figure 1, is a percolation-like transition above which a giant cluster of occupied sites percolates in space [37, 38]. More specifically, this second transition is—as it will be shown below—a standard/isotropic percolation phase transition 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 does not vanish [32, 33]. 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 (see figure 1). Let us stress once more, 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.

Hereafter, we study the nature of this (isotropic) percolation transition emerging within the active phase of the standard, or 'pure', contact process. However, our main focus will be on a variant of the contact process implementing external/environmental variability in the model, which is defined as follows: in the temporally-disordered contact process, at each Monte Carlo timestep the reproduction probability, p, takes a uniformly-distributed random value in the interval $\left[\overline{p}+\sigma ,\overline{p}-\sigma \right]$; thus, the parameter p of the pure contact process becomes a time-dependent function, p(t).⁶ The mean value, $\overline{p}$ is the control parameter and is kept constant in each simulation. In all cases, $\overline{p}$ is constrained to lie within the interval [σ, 1 − σ] so that p(t) ∈ [0, 1]; on the other hand, σ, characterizing the environmental variability, is typically fixed (unless explicitly stated) to a value σ = 0.25.

We start our analyses by focusing on the percolation transition of the 'pure' contact process—without additional heterogeneity—as it will be used as a benchmark for further analyses. The fundamental difference from the standard/isotropic percolation model is that the contact process and its variants are dynamical models and the states of nearby sites are correlated. Thus, the question we address now is whether such a difference is a relevant one, i.e. whether the universal aspects of the percolation transition observed within the active phase of the contact process are equivalent to those of isotropic percolation despite the existence of site-to-site correlations in the first case. In order to provide an answer to this question we resort to both, computational analyses and analytical/heuristic arguments.

The first observation to be made is that in the pure model there is a one-to-one correspondence between the control parameter p and the observed averaged density ρ (see figure 1 black curve), and as a result they are interchangeable⁷. Thus, we rely on the average density ρ in the forthcoming analyses, since it bears closer resemblance with what empirically measured in field analyses.

Figures 2(A) and (B) show the cluster-size distribution P(s) for a wide range of such occupation densities ρ in the pure contact process as measured in extensive simulations. At the percolation threshold p ≈ 0.725—corresponding to the critical occupation-density value ρ_c ≈ 0.56—and only at this value, the cluster-size distribution decays as a power-law, i.e. clusters of broadly different sizes emerge. Furthermore, the scaling behavior is fully compatible with the expected universal behavior of standard/isotropic percolation; indeed, the fitted value for the cluster-size distribution exponent τ in the contact process at the percolation threshold is in excellent agreement with the above-mentioned value τ ≈ 2.05 of standard isotropic percolation [32, 33]. Moreover, 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

$$\begin{equation}\chi =\sum _{i}{s}_{i}^{2}/{n}_{\mathrm{o}\mathrm{c}\mathrm{c}},\end{equation} \tag{ 4 }$$

—where i runs over all clusters, s_i is the corresponding cluster size and n_occ is the total number of occupied sites—is expected to obey the general scaling law

$$\begin{equation}\chi ={{\Delta}}^{-\gamma }\tilde {\mathcal{F}}\left({\Delta}{L}^{1/\nu }\right),\end{equation} \tag{ 5 }$$

where ${\Delta}=\vert {p}^{\mathrm{o}\mathrm{c}\mathrm{c}}-{p}_{\text{c}}^{\mathrm{o}\mathrm{c}\mathrm{c}}\vert$ is the distance to the critical occupation probability ${p}_{\text{c}}^{\mathrm{o}\mathrm{c}\mathrm{c}}$, L the system linear size, and ν = 4/3 and γ = 43/18 are exactly-known critical exponents for the two-dimensional isotropic percolation class [32, 33]. This implies that there exists a scaling function, o master curve, $\tilde {\mathcal{F}}\left({\Delta}{L}^{1/\nu }\right)$ to which the data should collapse for different system sizes L and Δ values (if the right values are used for the exponents ν and γ). Figure 2(C) shows such a curve collapse for the pure contact process at different values of p—around the critical percolation point p ≈ 0.725 (corresponding to ρ_c ≈ 0.56). Together, the three panels of figure 2 show that the pure/standard contact process (i.e. with σ = 0) displays a scaling behavior for a critical density ρ_c ≈ 0.56 that is fully compatible with that of the standard/isotropic percolation class.

Zoom In Zoom Out Reset image size

In order to rationalize this conclusion from a theoretical viewpoint, one can argue as follows. At the percolation critical point—which lies well inside the active of the pure contact process—occupied sites are not randomly placed as they are generated by a dynamical process that generates spatial correlations [31]. However, as the dynamics at such a density is far away from the absorbing-active phase transition (where long-range correlations are well-known to emerge), spatial correlations are only short-ranged. As a consequence—employing a renormalization-group perspective—at sufficiently 'coarse-grained' scales, the effective density becomes uncorrelated from one region to another. Thus, standard/isotropic percolation scaling needs to emerge at sufficiently large scales. In other words, the presence of weak short-range correlations is expected to affect non-universal features such as the location of the critical occupation density, but does not affect the large-scale universal behavior of percolation, in perfect agreement with our computational observations.

Summing up, the pure contact process 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 decay (with a characteristic scale) is observed.

When temporal heterogeneity (σ > 0) is introduced in the model, the scenario described in the previous section for the case of the pure contact process (σ = 0) is significantly altered, as we show in what follows. For the sake of simplicity and without loss of generality, we present results only for σ = 0.25, but analyses with other values have also been performed to verify the robustness of our conclusions.

The first important observation to be made is that—opposite to the pure version of the model—now there is not a one-to-one correspondence between the control parameter $\overline{p}$ and the resulting density ρ. Actually, as illustrated in figure 3, fixing $\overline{p}$ one observes a probability distribution of ρ values that does not become sharp, even in the infinitely-large system-size limit⁸. Owing to the intrinsic time variability of the control parameter, p(t), the systems keeps on jumping between different values of p—each one with a preferred stationary density—and, thus, the overall density fluctuates—without relaxing to a precise fixed value—even in the thermodynamic limit. This is also reflected in figure 3(A) in which we depict the resulting averaged density for each $\overline{p}$ (orange curve) together with its standard deviation (shaded region).

Zoom In Zoom Out Reset image size

As described e.g. in [49], the contact process with temporal disorder exhibits also an absorbing-to-active phase transition. Here, we find computationally that, for σ = 0.25, the absorbing-to-active critical point is located at ${\overline{p}}_{\text{c}}\approx 0.6550$, a value slightly larger than the corresponding value for the homogeneous (or 'pure') case p_c ≈ 0.6225, as illustrated in figure 3); for larger values of σ one can observe progressively larger shifts in the location of the critical point [49]. Let us also mention, that this temporally-disordered model exhibits an unusual scaling behavior, including logarithmic behavior for some quantities (see [49, 50]).

As already explained, as soon as σ > 0 the distribution of possible values of the overall density ρ is not a delta-function around a given density value (see figure 3(B)). A straightforward consequence of this is that the 'constant-$\overline{p}$ ensemble' and 'constant-ρ ensemble' are not equivalent: it is not the same thing to collect statistics of cluster sizes for a fixed value of $\overline{p}$ than to do so for a fixed value of ρ. This is in stark contrast with the standard pure case in which both ensembles are equivalent, as in the thermodynamic limit the distribution of ρ possible values becomes a delta function for any value of $\overline{p}$ [31, 49, 51].

In order to properly compare the computational results of the temporally-disordered model with field observations of vegetation cover—for which the only empirically measured quantity is the corresponding vegetation density ρ [17, 18]—it is more pertinent to consider a 'constant-ρ' ensemble. In practical terms, to obtain results for a given density ρ, we run computer simulations for different (homogeneously sampled) values of $\overline{p}$ and for sufficiently long times (typically t > 10³) as to reach a statistical steady state. Using such timeseries, we take 'snapshots' of different configurations across time (taking them sufficiently separated in time as to be uncorrelated) and collect the statistics of patches sizes conditioned to the value of the overall density ρ (±0.0005). Notice also that if an infinite (system-spaning) cluster exists, it is not considered for the statistics of finite clusters [33].

In what follows, we present the results of computational analyses for the resulting P(s) distributions for different values of ρ. Let us underline before going forward that, this way of proceeding considers configurations ('snapshots') generated under different dynamical conditions to generate one single probability distribution function for each density value. In other words, the obtained results for a fixed value of ρ correspond to a mixture of different dynamical processes, e.g. with different values of $\overline{p}$ and with diverse realizations of the temporal disorder, but nevertheless, well-defined probability-distribution functions emerge.

Using the statistics collected this way, one can find that the best fit to P(s) for each value of ρ. Given the subtleties usually encountered when determining whether a given distribution follows a power law, let us stress that in the forthcoming analyses we employ the standard most-stringent statistical tests (following the methods introduced in [52]; see appendix) to fit power laws and to decide whether they constitute a better fit than e.g. exponential distributions. Careful analyses reveal that—for a large variety of values of ρ—cluster-size distributions conform to power-law curves for at least two (and up to four) decades in size (see figures 4(A) and (B), and table 1). In particular, one can estimate the parameters C, τ, s_min, and s_max in equation (1) that maximize the normalized log-likelihood of the fit (see appendix).

Zoom In Zoom Out Reset image size

We observe the presence of fat tails in the distributions P(s) for different densities (observe also, that there are crossovers, clearly visible e.g. for ρ = 0.58). The tails of these distributions can be fitted as power-law functions with continuously varying exponents for a wide range of ρ values. Indeed, eye inspection of figures 4(A) and (B) already suggests that the exponent τ, characterizing the power-law decay, displays a non-monotonic variation with population density ρ (see also figure 4(C)). The fitted value decreases from τ ≃ 3.02 (for ρ = 0.02) down to τ ≃ 2.01 (for ρ ≈ 0.58) and increases again up to τ ≃ 3.53 (for ρ ≈ 0.7) (see figure 4(C) and table 1).

Moreover, in order to verify the robustness of these results, in figure 5 we present results for increasing system sizes, confirming that the measured power-law fits are not finite-size dependent, i.e. they are robust. Notice, in all cases, the lack of strong downward bends even if the considered system sizes are not very large; i.e. the distributions are consistently fat tailed.

Zoom In Zoom Out Reset image size

Thus, remarkably, opposed to the case of the pure contact process—power-laws fit well the cluster-size distributions P(s) not only at a unique (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, at least, up to numerical resolution. Owing to this peculiarity, the percolation threshold ρ_c cannot be easily obtained—as usually done in numerical analyses—just by looking for the precise value of the control parameter for which scale-free cluster-size distributions emerge [49]. Therefore, in order to obtain an estimate of the critical point it is necessary to employ some criterion; in particular, we define the critical percolation point as given by the value of ρ for which τ changes its tendency, from decreasing to increasing, i.e. ρ ≈ 0.58 (see table 1 for further details on the fits). We need to emphasize that this is only an approximation; more precise estimations of the percolation critical point are not easy to obtain, not even employing finite-size scaling analyses [49, 50].

The cluster-size distribution exponent fitted at ρ = 0.58 is τ ≃ 2.01, slightly smaller than the expected value in two-dimensional percolation τ ≈ 2.05 (see table 1). Moreover, for the present estimation of the critical point, we obtain a scaling collapse (see figure 4(D)) similar to that of the homogeneous case, represented in figure 2(C) but with exponents ν ≃ 3 and γ ≃ 4 that are considerably larger than the percolation values (ν = 4/3 ≈ 1.33... and γ = 43/18 ≈ 2.39...), suggesting a different type of critical behavior (or universality class). However, it is important to remark that by slightly shifting the critical point estimation, e.g. to ρ = 0.56, for which τ ≈ 2.31, we obtain a similar-quality curve collapse with slightly different exponents (ν ≃ 3.5 and γ ≃ 4). This suggests the existence of relatively large error bars in the numerically-estimated critical exponents.

Thus, at the light of all these difficulties, it is clear that a more precise and careful determination of the critical point, the associated critical exponents, and their possible scaling relationships is still needed to obtain a deeper understanding of this novel type of percolation phase transition. This is, however, beyond the scope of the present manuscript and is left as a challenge for future work.

We have shown that the contact process equipped with temporal disorder exhibits scale-free cluster-size distributions without the need of careful fine tuning to a threshold point. Nevertheless, in order to be able to make a more direct comparison with empirical observations, an important caveat remains. The results obtained so far correspond to ensemble averages, where many different realizations of the system—e.g. at different times—are considered together with the only restriction of a fixed ρ value. However, empirically measured exponents of vegetation cluster-size distributions are determined for individual 'snapshots'.

To bridge this gap, we 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 (homogeneously sampled) values of $\overline{p}$ and took instantaneous pictures of the system configuration at the steady state. For each of such snapshot, j, we computed its overall density and measured its associated cluster-size distribution (histogram), P_j(s). Obviously, these distributions—not being ensemble averaged—are much noisier than their ensemble-averaged counterparts. However, simple eye inspection of figure 6(B) already reveals that the cluster-size distribution of an instance, for an individual snapshot, can span over a large range of cluster sizes, and can, in fact, be fitted by a power-law, even for a density far away from the percolation threshold (such as ρ = 0.3 in this example).

Zoom In Zoom Out Reset image size

A collection of individual-snapshot cluster-size distributions are shown in figure 7 for two different densities: ρ = 0.48 and ρ = 0.58. This figure illustrates that even for a fixed density ρ there is a high variability in the shape of 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. Remarkably, even if exponentially decaying curves do exist, especially for low densities, power-laws are present in snapshots for most densities. These power-law fits are more prominent and extend for a broader size-interval in the neighborhood of the percolation transition (see table 2 for more details on the fits).

Zoom In Zoom Out Reset image size

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 $\overline{\tau }\left(\rho \right)$; 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 (i.e. power-laws give a higher normalized log-likelihood). The fifth column shows the fraction for which the power-law fit extends for at least 1.5 decades.

Density	$\overline{\tau }$	Exponent	% power law better	% more than 1.5
	exponent	variance	than exponential	decades fit
0.04	2.14	1.10	71.74	71.74
0.06	2.38	1.26	82.54	61.90
0.08	2.46	1.21	81.48	59.26
0.10	2.63	1.33	78.95	40.00
0.12	2.67	1.42	70.16	32.26
0.14	2.75	1.32	64.29	26.62
0.16	2.79	1.30	62.23	31.38
0.18	2.70	1.34	62.00	28.50
0.20	2.57	1.24	73.50	32.50
0.22	2.49	1.24	81.00	43.00
0.24	2.35	1.13	93.00	37.00
0.26	2.43	1.36	100.00	45.50
0.28	2.33	1.24	100.00	55.00
0.30	2.16	1.05	100.00	54.00
0.32	2.05	0.89	100.00	60.00
0.34	1.93	0.63	100.00	66.50
0.36	1.80	0.37	100.00	70.00
0.38	1.71	0.23	100.00	75.50
0.40	1.62	0.11	100.00	84.00
0.42	1.58	0.07	100.00	87.50
0.44	1.60	0.07	100.00	92.00
0.46	1.62	0.06	100.00	94.97
0.48	1.66	0.07	100.00	94.95
0.50	1.70	0.06	100.00	96.97
0.52	1.77	0.07	100.00	96.48
0.54	1.86	0.11	99.49	95.43
0.56	1.96	0.18	100.00	93.91
0.58	2.01	0.15	100.00	92.89
0.60	2.10	0.18	100.00	89.23
0.62	2.15	0.23	99.48	88.02
0.64	2.24	0.32	98.97	63.92
0.66	2.39	0.57	99.47	61.38
0.68	2.50	0.71	97.30	80.54

We also compared the quality of fits to power-law functions with those to e.g. exponential-distribution P(s) = Ce^−αs. For many snapshots—especially when close to the percolation critical point—power-law distributions constitute a better fit to the data than exponential ones (i.e. power-laws give a higher normalized log-likelihood) but this is not always the case, especially for low densities ρ ∈ [0.04, 0.24] (see fourth column of table 2 and appendix for further details). Nevertheless, even for densities in this interval, power-law fits are better in at least the 60% of the cases (see table 2 and appendix). On the other hand, for intermediate densities, power-law fits outperform exponential ones in almost each one of the snapshots. Observe that in spite of the existing limitations (i.e. the finite system size and the lack of an averaging procedure), power-laws are observed for at least 1.5 decades in logarithmic scale in most of the cases (fifth column of table 2).

To avoid confusion, let us underline that some of these effective exponents are smaller than 2 and that power laws with exponents smaller than 2 do not have a well-defined averaged value when integrated to arbitrarily large values of the variable. This means that the fits are just approximated ones and cannot possibly extend to arbitrarily large cluster sizes (see [53] for an illuminating discussion of this and related issues).

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 figure 8) similar to that of the average in figure 4(C).

Zoom In Zoom Out Reset image size

It is important to emphasize that this effective averaged value of $\overline{\tau }\left(\rho \right)$ is not the same as the ensemble averaged value of τ(ρ) obtained in the previous section: here, we fit individual snapshots, discard those that comply better with exponential fits, and average the resulting set of τ values to obtain a final averaged effective exponent $\overline{\tau }\left(\rho \right)$; this is different from considering the averaged distribution including all snapshots with a given ρ, and fitting a single exponent τ(ρ) to such a final distribution. Observe that the mean exponent $\overline{\tau }\left(\rho \right)$ reaches a minimum value $\overline{\tau }\simeq 1.58$ at a density of ρ ≃ 0.42, below the percolation threshold, and increases up to $\overline{\tau }\simeq 2.8$ 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; in the limit in which σ goes to 0 one obtains typically exponential fits, and one should recover power-law scaling only in a vanishing region around the critical point.

Finally, we can compare the obtained averaged exponents $\overline{\tau }$ with available empirical data [17, 18]. In particular, the black dots of figure 8 correspond to empirical data for different vegetation-cover densities. Remarkably, field data are roughly compatible with the model results (for σ = 0.25) within one standard deviation (p-value p = 0.04 < 0.05 of Pearson's χ² statistics). Of course, this can be just a coincidence, as we have no information on what the empirical value of σ would be for the different actual data, nor it needs to be same for all of them. Still, we find it noteworthy that the overall dependency on vegetation density is reproduced by the present simple model.

The study of vegetation patterns has a long tradition in ecology and environmental sciences. This interest has been reinforced in recent years owing to the possibility 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 [17, 18]. Puzzled by these field observations, theoreticians were challenged to understand the origin of such scale-free patterns and, as a result, diverse mechanisms have been proposed to account for them. Different models put the emphasis on existing plant–plant interactions (competition and facilitation), feedback loops between vegetation dynamics and abiotic resources such as water [17, 19], or the presence of spatial and/or temporal heterogeneity in the environment [18, 20, 21, 54].

For the sake of generality, let us stress here that the theoretical challenge of modeling 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 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 conundrum 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 [55–58], 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 [59]. Such type of models already exist in statistical mechanics (for an in-depth discussion of these issues we refer to [59]), 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 [35, 36, 60–62]. 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 the emergence of very broadly distributed vegetation patterns—whose tails can be fitted as power laws, i.e. the hallmark of scale invariance—for a very wide range of parameter values. Our model is individual-based and spatially explicit, based on the two-dimensional contact process, an extensively studied model that is the simplest to mimic vegetation dynamics (as a birth-and-death process). The contact process 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 average density. In other words, the contact process 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 figures 1 and 2) 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 [37, 38] 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 (see figure 2). 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 renormalization-group theory could have correctly anticipated though—employing renormalization-group reasonings—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 [49, 50, 63], 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 justified—novel ingredient the model changes dramatically.

• (a)  
  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), even if further analyses are required to establish the precise values of the exponents and the associated scaling relationships.
• (b)  
  Second, generic scale-free distributed clusters emerge on average for a broad range of overall densities and not just at the critical one. As in the previous point, further work is needed to precisely establish the values of their asymptotic exponent values and their dependence on ρ.
• (c)  
  Third, individual snaphots of the system state—i.e. single configurations in which no averaging is performed—can also exhibit clusters of highly variable sizes that can be fitted as power laws for a broad range of possible vegetation densities. Even more remarkably, the obtained power-law exponents are density dependent—in agreement with what observed in empirical findings—and exhibit—roughly speaking—the same type of density dependence as those observed empirically.

Let us now briefly discuss each of these three results.

• (a)  
  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, 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 beyond our scope here. However, let us remark, that usually the presence of temporal disorder dramatically affects the universal behavior of well known universality classes [48, 50] (including the quiescent-active critical point of the contact process, which is dramatically affected by temporal noise [49, 64, 65]). 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?) type of scaling, establishing all critical exponents and scaling relationships with good accuracy and precision—and determining whether this really constitutes a robust universality class—is left for future work.
• (b)  
  The emergence of a broad region of densities for which scale-free behavior 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 [49] in analogy with the well-known Griffiths phases that emerge in the presence of spatial (quenched) heterogeneity [66]. In particular, it was shown in [49] (see also [63, 72]) that for models similar to the contact process studied here, some 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, and exhibit an 'exotic' type of scaling [50]).
• (c)  
  Directly related to our findings is the important work of Manor and Shnerb who discussed—in a mean-field setup different from our spatially explicit one—how multiplicative (environmental) noise can lead to scale-free distributed cluster sizes [19, 21]. 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) [27, 67, 68]. 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, we believe that the reason why our model produces generic power-law distributions (at least in an approximate way) is twofold: (i) on the one hand, nearby the standard percolation phase transition clusters of rather different sizes are generated and (ii) owing to the multiplicative process induced by changing environmental conditions such clusters 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 [69]. In any case, further work is still needed to confirm in a more clean cut way that the power laws observed here within a broad range of parameter values are indeed asymptotic 'true' ones, observable at arbitrarily large scales; deeper theoretical understanding is still needed to advance in this direction.

The conclusion of the present study is that a very simple individual-based model with fluctuating external conditions is able to reproduce 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, in agreement with what empirically observed.

From the theoretical-ecology perspective, 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 confirms that temporal variability is a key player to explain the emergence of scale-free vegetation patterns in semiarid environments.

From a theoretical-physics perspective, a new type of percolation transition has been uncovered here. It is our hope that this works stimulates further analyses of its scaling properties and helps to develop new applications to real-world phenomena.

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 ERDF, Ref. A-FQM-175-UGR18 and SOMM17/6105/UGR and for financial support. We thank Juan A Bonachela, Paolo Moretti, Victor Buendía, Guillermo Barrios, Johannnes Zierenberg, and Pablo Villegas for a critical reading of the manuscript as well as for very useful comments.

We fit cluster-size distributions P(s)—both averaged over different realizations and times and for individual snapshots—as power-laws P(s) = Cs^τ within an optimal interval [s_min, s_max]. For each combination of the set of parameters τ, C, s_min, s_max, we obtained 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 [52, 70, 71]. We then selected the parameter sets that maximizes the normalized log-likelihood in each case. To avoid overfitting, we discard fits whose optimal range [s_min, s_max] includes less than a fixed number of points: for averaged distributions P(s) (figure 4) this range is fixed to 25 points, while for snapshots (figure 7) we considered intervals of more than 10 points. Also, to ensure a fitting interval [s_min, s_max] within the asymptotic region [53] of averaged distributions, we limited s_min to be large enough as to overcome crossovers at small-size clusters. In all cases, we compared the quality of power-law fits with the corresponding one for exponential ones, P(s) = Ce^−cs, within the obtained interval [s_min, s_max].