Quantitative analysis of self-organized patterns in ombrotrophic peatlands.

We numerically investigate a diffusion-reaction model of an ombrotrophic peatland implementing a Turing instability relying on nutrient accumulation. We propose a systematic and quantitative sorting of the vegetation patterns, based on the statistical analysis of the numbers and filling factor of clusters of both Sphagnum mosses and vascular plants. In particular, we define the transition from Sphagnum-percolating to vascular plant-percolating patterns as the nutrient availability is increased. Our pattern sorting allows us to characterize the peatland pattern stability under climate stress, including strong drought.

Introduction

Peatlands are humid ecosystems, poor in O2 and often in nutrients, with over 65% of organic contents over dry soil mass[1]. Covering around 4000000 km2 worldwide, they constitute an essential carbon well storing 270 to 550 PgC[2], an amount comparable to that of carbon released into the atmosphere since 1750[3]. Besides their contribution to the carbon cycle, the stratified accumulation of organic matter in a relatively preserved environment constitutes remarkable archives of past climate. Their dynamics is mainly governed by their hydrology[1,4]. Ombrotrophic peatlands, that are isolated from surface and subsurface water flows, primary rely on precipitations[5].

First attempts to model peatlands primarily aimed at describing peat formation, with a focus on the vertical dimension[6-9]. However, like many ecosystems[10,11], peatlands display self-organized spatial patterning, that is generally considered stable although it may slowly evolve in time[12]. The peatland patterns are believed to be resilient against climate changes[4,13-15]. Based on two-dimensional modelling, several processes have been proposed to create or influence these patterns, including peat accumulation, water ponding, and nutrient accumulation[9,16,17]. In the following, we shall focus on the latter, based on numerical[16] and empirical[18] arguments. Nutrient accumulation can be described as a reaction-diffusion process, in a two-dimensional model. The active transport of nutrients towards spots of vascular plant like Carex, results in their accumulation at specific locations where the subsequent development of vascular plants is favored, and in their depletion elsewhere[19]. This model was refined by considering the competition between vascular plants and Sphagnum mosses[20], leading to sharper and more stable patterns.

Here, we analyze the biomass-nutrient interaction in terms of a Turing instability[21-23], relying on the strong discrepancy between the fast and slow diffusion of nutrients and vascular plant, respectively. We further build on this class of two-dimensional reaction-diffusion model to systematically analyze and classify the resulting patterns, based on statistics of the number and filling factor of clusters of vascular plants and Sphagnum. We also investigate the equivalent phase transition from Sphagnum and vascular plant cluster percolation, and characterize the associated scaling law. Finally, we apply our classification to describe the peatland pattern stability especially under climate stress, including strong prolonged drought.

Peatland Model Description

Our work is based on Eppinga et al.’s two-dimensional reaction-diffusion model[20], which describes the interactions between vascular plants, Sphagnum, nutrient availability and hydraulic head (Fig. 1). The heart of the interaction is nutrient accumulation due to water draining by vascular plants, via transpiration. Vascular plants also promote water loss by their transpiration and compete for space with Sphagnum.

Figure 1

Main interactions between the four sub-systems of the peatland model.

More specifically, the model relies on four partial differential equations describing the temporal evolution of the following spatial variables: vascular plant and Sphagnum mosses biomass per unit surface and , hydraulic head , and nutrient availability , the latter being linked to the nutrient concentration per unit surface (both in g/m2) bywhere θ describes the soil porosity. The vascular plant biomass density evolves as:where the first term is the growth of the plants, rB being a growth parameter and the water-stress function f(H) defines water availability (f(H) = 0 if the water level is below the vascular plant root depth, 1 if above a threshold level h1 + z, and evolves linearly in-between). The second term encodes the vascular losses through death (d) and other causes (b). The third and fourth terms represent local environmental modifications due to the presence of Sphagnum mosses, as detailed in[20]. The last term accounts for diffusion. Table 1 summarises the meaning and default value of each variable.

Table 1

Parameter and variables symbols, units and default values[10,20].

Symbol	Interpretation	Unit	Default value
Variables
B	Vascular plant piomass	gB m−2	Variable
S	Sphagnum biomass	gS m−2	Variable
H	Hydraulic head	m	Variable
N	Nutrient availability	gN m−2	Variable
Nc	Nutrient concentration	gN m−3	Variable
water-stress function: f(H)=H−z−h_2h_1−h_2 for h1 ≥ H − z ≥ h2, f(H ≤ h2 + z) = 0, f(H ≥ h1 + z) = 1
h	Pressure head	m	Variable
h1	Pressure head below which soil water stress occurs	m	0
h2	Root depth	m	−0.5
z	Elevation head	m	1
Vascular biomass equation
rB	Vascular plant growth parameter	m3gN−1 yr−1	0.2
d	Vascular plant recycling parameter	yr−1	0.1
b	Vascular plant loss parameter	yr−1	0.2
DB	Vascular plant diffusion parameter	m2 yr−1	2
cBS	Sphagnum inhibiting coefficient on vascular biomass	gBgS−1yr−1	0.1
kBS	Feedback coefficient for vascular plants	gBgS−1 yr−1	0.5
HB	Half-saturation density for vascular plants	gB m−2	800
Sphagnum biomass equation
rS	Sphagnum growth parameter	yr−1	0.2
Smax	Maximum density of Sphagnum mosses	gSm−2	800
kSB	Feedback coefficient for Sphagnum	gSgB−1 yr−1	0.1
HS	Half-saturation density for vascular plants	gS m−2	300
DS	Diffusion coefficient for Sphagnum	m2 yr−1	0.2
Hydraulic head equation
p	Precipitation	m yr−1	0.5*
tv	Transpiration parameter	m3gB−1 yr−1	0.005*
e	Evaporation parameter	m yr−1	0.3*
k	Hydraulic conductivity	m yr−1	500
θ	Soil porosity	Dimensionless	0.7
Nutrient availability equation
Nin	Nutrient input	gN m−2 yr−1	0–5
u	Vascular plant uptake parameter	m3gB−1 yr−1	0.002
r	Nutrient loss parameter	yr−1	0.1
DN	Diffusion coefficient for nutrients	m2 yr−1	10

*Denotes values affected by climate. gB, gS, and gN respectively denote grams of vascular biomass, Sphagnum biomass, and nutrients.

The growth of vascular plants is regulated by their access to nutrients. In an ombrotrophic peatland, the nutrient external input Nin mostly originates from the atmosphere, via precipitation[1,24]. The nutrient availability N is then governed by losses through the uptake by the vascular plants (with u the uptake parameter), the release of nutrients by dying vascular biomass, other losses characterised by the parameter r, and nutrient transport allowing the vascular plants to harvest nutrients from their neighbourhood. This transport occurs both by diffusion according to Fick’s law and through advection with the groundwater. The latter is described by the Darcy flow, which expresses fluid flows at Reynolds numbers below 1–10, i.e. the regime where inertial effects are negligible, so that the flow occurs along the pressure gradient. This regime is typical of geophysical flows, including the diffusion of water through porous materials. Thus, the nutrient availability evolves aswhere k is the hydraulic conductivity. The vascular plant access to nutrients is modulated by the access to water. The hydraulic head H evolves under the influence of the precipitations p (with soil porosity parameter θ), the evaporation parameter e, and the transpiration parameter tv:

Note that the transpiration only depends on the vascular plant biomass B, as Sphagnum mosses have no vascular system and do not transpire. The last term is the Darcy Flow. Finally, vascular plants compete for space with Sphagnum mosses, the biomass S of which is given bywhere the first term describes a logistic growth, while the second one corresponds to competition with vascular plants, and the last one is diffusion. All parameters driving the model are spatially homogeneous and constant. They are defined as mean annual values, without seasonal nor daily variability. The numerical implementation is detailed in the Methods section below.

It should also be pointed out that the model primarily aims at characterizing the self-organized patterns, but does not seek perfect realism. In particular, the model is restricted to two types of plants, namely mosses and a generic vascular plant, not distinguishing between, e.g., trees, shrubs, and sedges[1]. The evolution of the soil, including peat formation, decomposition, and loss[25] and the associated deviation of the peatland surface from a planar one is also disregarded, preventing the model from realistically describing the differential soil porosity and conductivity to water[26-28]. Finally, the periodic boundary conditions and the flat topography exclude any net water flow through and into the simulation area[5], and restrict water input to precipitations[17].

The effect of climatic stress was investigated by plugging the data from the 2 × CO2 scenario of the Canadian Climate Center Second generation General Circulation Model (GCMII)[29]. We selected the simulation results from a peatland-rich area in Siberia (grid point of 3.75° × 3.75° size centred on 90° East, 61.23° North) and we used the yearly averages of the precipitation, as well as the evaporation and transpiration parameters from the GCM as inputs for our simulations. Conversely, the dependence with the climate scenarios of the nutrient input rate Nin was not considered in the present work.

Results and Discussion

Pattern analysis

Figure 2 illustrates the variety of peatland patterns in the framework of the nutrient concentration-limited model, as already pointed out by Rietkerk et al.[19]. It displays the patterns observed after 400 years, once a stable regime is reached. Initial conditions are described in the Methods below and parameters in Table 1. Six qualitatively distinct patterns can be distinguished, for increasing external nutrient input Nin: (1) full Sphagnum coverage; (2) sparse quasi-circular spots of vascular plants; (3) lines and spots of vascular plants, interleaved with a percolating Sphagnum network; (4) percolating network of vascular plants; (5) dominance of vascular plants with isolated Sphagnum spots; and (6) full vascular plant coverage. Note that, in this work, we define percolation as the existence of at least one cluster continuously connecting two opposite sides of the simulation area. Although patterns form much earlier in the simulations (See Supplementary Movie 1), they keep evolving for long times and only stabilise after about 400years. Although we do not explore it here, the evolution of peatland over shorter time duration is also of interest (see for example Morris et al.[17] for evolution over a decadal range).

Figure 2

Simulated peatland patterns after 400 years stabilization. Initial conditions and parameters are as specified in the Methods section and in Table 1. From top to bottom: nutrient availability distribution, water level, Sphagnum biomass, and vascular plant biomass. Columns 1–6 correspond to nutrient inputs N = 0.25, 1.25, 2.25, 3, 3.875, and 4.5 g/m2/yr, respectively, as specified above each column. Each panel represents a square 256-meters wide area, with 2 m resolution. See Supplementary Movie 1 for the pattern build-up from random initial conditions for Nin = 2.5 g/m2/yr.

Comparing the patterns for Sphagnum and Vascular plant biomasses illustrates the mutual exclusion between the two species, which also explains the sharpness of the patterns. The anti-correlation between the vascular plant biomass and the hydraulic head is due to the water depletion by the vascular plant evapotranspiration. The nutrient availability pattern is slightly more complex. Nutrient maxima at the vascular plant hot spots are surrounded by a nutrient-depleted region, as is especially visible in the second panel of the first row. Due to this structure, vascular plant regions are self-attracting at short distance and self-repulsive at medium range, over a distance governed by the nutrient harvesting efficiency of the vascular plants. This distance governs both the size of the vascular plant clusters (either diameter, or line width), and the distance between them, i.e., the system quasi-periodicity. Furthermore, Regimes 1 and 6, 2 and 5, and 3 and 4, respectively, correspond to each other by exchanging the roles of vascular plants and Sphagnum. The variety of the patterns is typical of a Turing instability driven by a reaction-diffusion model with mutual retroactions and differential diffusion speeds[21-23]. Here, the fast diffusion of water and nutrients as compared to that of biomass provides the adequate contrast for generating dots, lines, branched structures and percolating clusters from the same species, depending on a single control parameter, Nin.

Beyond the qualitative description of the patterns and the understanding of the interactions between the model constituents, our aim is to provide a quantitative characterization and systematic sorting of the clusters of Sphagnum and of vascular plants, respectively. To this end, we plotted the number of clusters of both vascular plants and Sphagnum as a function of Nin (Fig. 3). In the low-nutrient regime (Regime 1) neither vascular plant clusters, nor non-percolating Sphagnum clusters are observed. The jump into Regime 2 corresponds to the point where nutrients are sufficient to allow the long-term survival of sparse vascular plant clusters, the number of which increases with the nutrient input Nin. Simultaneously, their biomass density increases and their nearest-neighbour distance decreases, allowing interactions between neighbouring clusters. They start merging into linear clusters, defining a second threshold above which the number of vascular plant clusters starts decreasing. Progressively, these linear clusters branch and tend to isolate Sphagnum regions, allowing the appearance of non-percolating Sphagnum clusters in increasing numbers.

Figure 3

(a) Number of clusters as a function of nutrient input. Insets: typical patterns of vascular plant biomass in each range, for N = 0.25, 1.25, 2.25, 3, 3.875, and 4.5 g/m2/yr. The vertical dotted lines at N = 0.31, 1.58, 2.63, 3.35, and 4.41 g/m2/yr correspond to the transitions between regimes. (b) Filling factor (with a binarization threshold at 500 and 200 g/m2 for vascular plant and Sphagnum biomass, respectively) as a function of the nutrient input. See Table 1 and Methods for parameter values and initial conditions.

The transition from Sphagnum percolation to vascular plant percolation defines the transition from Regimes 3 to 4. Although the percolation threshold is characterised independently from the cluster number (See Methods), it approximately corresponds to the crossing between the decaying number of vascular plant clusters, and the rising number of Sphagnum clusters.

The initially disconnected vascular plant clusters then progressively merge into the percolating one. Simultaneously, the number of Sphagnum clusters increases, as the Sphagnum domains divide into smaller entities in reaction to the emerging connexions between the vascular plant clusters. This number reaches a maximum where the Sphagnum clusters have their minimum stable size. In Regime 5 this number decreases again because the vascular plants diffuse faster and spread at the expense of Sphagnum clusters. When the latter are below their minimum stable size, they finally vanish. Note that during this process, the Sphagnum biomass density per unit surface where it is present is stable and close to its maximum value Smax, so that the Sphagnum decay mainly occurs via the reduction of the surface it covers. Finally, in Regime 6, Sphagnum has fully disappeared under the pressure of the vascular plants.

We characterized the percolation transition by determining the scaling law and the associated critical exponent near to it. We fitted the mean cluster area (excluding the percolating cluster) of both vascular plant and Sphagnum, as a power law of the filling factor p, i.e., the fraction of the surface covered by the considered species (Fig. 3b) in the vicinity of the transition. This transition occurs at the critical filling factor pc (pc,V = 0.37 and pc,S = 0.63 for binarization thresholds of 500 and 200 g/m2 for vascular and Sphagnum biomass, respectively):

We obtain and for the vascular plant and Sphanum biomass, respectively (Fig. 4). Although the filling factors depend on the binarization threshold, the latter has little influence on the critical exponents within the ranges discussed in the Methods below. These values are much lower than observed in other experimental percolating systems[30] as well as the theoretical prediction of for two-dimensional uncorrelated percolation. This smoother transition to percolation in our model expectedly stems from local correlations due to the plant growth and expansion dynamics, as well as the various feedbacks between the biomass, hydraulic head and nutrient concentration.

Figure 4

Power scaling of the average cluster area in the vicinity of the percolation transition. The filling factor is calculated by binarizing the patterns at 500 and 200 g/m2 for vascular plant and Sphagnum biomass, respectively.

Effect of climate on peatland patterning

The above-described sorting of self-organized patterns and characterization of the transitions between them allows to investigate the pattern changes and therefore their stability more quantitatively. Here, we illustrate this ability in the case of climate change.

The model parameters that may be expected to be sensitive to climate change are the ones related to the hydrological fluxes: precipitations, evaporation, and transpiration. The latter two strongly depend on temperature. For example, the evapotranspiration is usually expressed by the Penman-Monteith equation[31]. Moreover, the evaporation and transpiration terms of Equation (4) display similar trends when atmospheric conditions change. Experimental studies showed that both are mainly governed by the water availability in the soil, i.e., the water-stress function f(H). In particular, the transpiration parameter tv is quite insensitive to changes in temperature: in hot conditions transpiration cooling tends to compensate the reduced evaporation related to photosynthesis[32,33]. Climate change is therefore implemented into the model by changing the precipitation level as well as by multiplying the evapotranspiration term by a Penman-Monteith factor that is directly provided by the GCMII.

We first performed a single-parameter sensitivity analysis to estimate the model response to a change in the climatic variables. As displayed in Fig. 5a,b, the transitions between the six regimes depend little on precipitation and evaporation, within a wide range of realistic values. These parameters therefore have little impact, if any, on the peatland patterns in the framework of the model. Beyond the visual observation of similar patterns, our approach is more quantitative and supports an interpretation in terms of nutrient accumulation of previous assumptions that peatlands are robust to climate change[13-15]. This stability also illustrates the fact that within our peatland model, changes in net water supply are not critical, while the parameter primarily governing the pattern formation in this system is nutrient availability and its local accumulation by vascular plants. In contrast, strongly decreasing the vascular plant transpiration parameter tends to restrict the range of nutrient input over which the transitions occur (Fig. 5c). This can be understood by considering that a reduced vascular plant transpiration will reduce the hydraulic head gradients and the associated nutrient transport via the groundwater flow (see last term in Equation (3)). At low external nutrient input Nin, the reduced nutrient flow towards vascular plants limits their ability to harvest nutrients. Therefore, the threshold for their survival increases. At higher nutrient input, vascular plants can develop, but the less efficient transport of nutrients limits the inhomogeneities in their concentration, hence also in the vascular plant distribution. The latter therefore spread more, so that the thresholds for their domination over Sphagnum is lower. However, the transpiration parameter tv depends mostly on the plant species and the water-stress function f(H), rather than on atmospheric temperature and humidity. In particular, Kim and Verma[32] showed little deviation between the measured and potential evaporation. Furthermore, Drake et al. observed in Eucalyptus that a simulated heat wave induces a simultaneous decrease of photosynthesis and increase of transpiration cooling, with a conserved total water flux, governed by the water availability[33]. Therefore, a significant change in evapotranspiration would require a change in the vascular plant species, which is beyond the scope of the present study.

Figure 5

Thresholds of the self-patterning regimes as a function of (a) precipitation, (b) evaporation parameter e, and (c) transpiration parameter tv. See Table 1 and Methods section for parameter values and initial conditions.

To further investigate the stability of peatland patterning under climate stress beyond the variation of individual parameters, we simulated a change in the climate. After stabilizing the patterns over 400 years, we switched the precipitation and evapotranspiration rates to follow those from the 2 × CO2 scenario of the GCMII. As expected from the stability of the thresholds against changes in individual parameters, the patterns remain very stable in response to the climate change (See the behavior before year 415 and after year 424, i.e., outside of the drought time, in Supplementary Movie 2). This is observed over a wide range of nutrient input, namely from 0.25 to 4.75 g m−2 yr−1. The system equilibrates by adapting the vascular plant biomass (Fig. 6a,b) to the available water, stabilizing the hydraulic head gradient and level by increasing the water losses by transpiration to compensate the precipitation level changes (Fig. 6c,d). This feedback is faster (at the year-scale) for high nutrient availability (), and slows down when the input nutrient flux is lower (see Fig. 6). There are two thresholds in the response. The first one corresponds to the resilience of vascular biomass (i.e., the transition from Regimes 1 to 2, see Figs 2 and 3). Without vascular biomass (Regime 1), evapotranspiration corresponds to evaporation alone and the total water flux is not balanced. From Regime 2 on, the system responds (Fig. 6a,c) with a certain time-lag. For , there is enough vascular biomass in the system to respond rapidly (yearly time scale) to variations in precipitation (Fig. 6b,d). This feedback ensures the resilience of the system. By influencing the density of vascular biomass, it will also expectedly influence the efficiency of carbon storage in the peatland.

Figure 6

Influence of the 2 × CO2 scenario of the GCMII on a patterned peatland. (a,b) Evolution of total vascular biomass on simulation area; (c,d) Precipitation and spatially-averaged evapotranspiration. Climate fluctuations are turned on after stabilising the patterns over 400 years. See also the corresponding pattern stability in Supplementary Movie 2.

The strength of these feedbacks even allows the patterns predicted by our model to promptly recover after an extreme drought. Figure 7 displays the effect on the spatial mean of vascular biomass and nutrients of a 10-years long drought reducing the precipitations from a mean value of 800 mm/year to 50 mm/year. The pattern structure is not affected by this severe drought (See years 415–424 in Supplementary Movie 2), but the level of biomass strongly decreases, with a drop by 67% after 10years. Correspondingly, the nutrient intake by vascular plants is reduced, so that nutrients accumulate in the soil during the drought (Fig. 7b). These extra nutrients allow an efficient recovery of the vascular plant biomass at the end of the drought, resulting in a 4-years-long overshoot of up to 35% above the steady-state vascular biomass density observed before the drought. Finally, once the stored nutrients are consumed, the system recovers to its initial state within 12 years of the end of the drought.

Figure 7

Peatland response to a 10-years long drought. (a) Temporal evolution of spatial mean of vascular biomass and precipitations; (b) Spatial mean of nutrient availability. The initial condition before the drought is a peatland stabilized after 400 years under the same conditions as in Fig. 2, with Nin = 1.375 g m−2 yr−1, i.e., in Regime 2. Climate conditions are then switched to the 2 × CO2 scenario of the GCMII, except for years 415–424 where a drought is simulated by reducing precipitation to 50 mm/year. See also Supplementary Movie 2 illustrating pattern resilience during the drought of years 415–424.

Similar behaviours are observed for various magnitudes and durations of the drought, as well as over a wide range of nutrient inputs (0.25–4.75 g m−2 yr−1, corresponding to the Regimes 2–6 displayed in Figs 2 and 3). The magnitude of the peatland response to a drought depends on the duration and the precipitation level of the drought, as well as on the nutrient input. On the other hand, the recovery time appears to be essentially independent from both the hydrology and nutrient input.

It should however be kept in mind that the model we use in this work is simplified. In particular, the hydrology is reduced to a homogeneous water table level. Hence, the impact of water shortage is handled via the water-stress function f in the evolution of the vascular plants, but the model does not allow for a full depletion of the water table and therefore does not describe the decay of the Sphagnum mosses biomass under dry conditions (see Equation (5)). Considering Sphagnum decay would expectedly result in a decrease of the pattern contrast[18], but without major influence on the rest of the model as the only coupling considered here is via the competition for space with vascular plants. Due to these limitations, the model outcome as described above cannot be expected to provide quantitative predictions, but rather indications about trends like the robustness of the nutrient accumulation mechanism, and a testbed for interpreting empirical observations. Biologically and ecologically realistic predictions would require a much more detailed and complex model, as discussed above.

Conclusion

As a conclusion, based on a diffusion-reaction model of an ombrotrophic peatland, we have proposed a detailed analysis of vegetation pattern, with a quantitative description of the pattern transitions and a systematic sorting of the self-organized patterns. We characterized the transition from Sphagnum to vascular plant cluster percolation when the nutrient availability is increased, as well as the associated scaling laws close to the transition. Such classification and characterization is not specific to the model of Eppinga et al.[20], but can straightforwardly be extended to any system displaying Turing instabilities, e.g., reaction-diffusion models with differential spatial mobilities. This covers many ecosystems[10,22], including non-ombrotrophic peatlands[34], mussel banks[11], coral reefs[35], drylands[36], or savanna[37].

Applying our pattern classification to the self-organized patterns generated by a peatland model based on nutrient accumulation to a situation of climate stress, including strong prolonged drought, we evidenced its high resilience. This finding may help understand the empirically observed high peatland stability and resilience[4,13-15].

Methods

Numerical implementation

Simulations were performed in Matlab on a 128 × 128 pixel grid with periodic (toroidal) boundary conditions, covering 256 × 256 m, using a Newton (Euler) numerical scheme with a temporal step dt = 10−4 year. We started from an initial random distribution of spots bearing 200 g m−2 biomass of either vascular plants or Sphagnum mosses, each covering 12.5% of the grid elements. We checked that the results are independent from the particular initial random pattern. The initial hydraulic head and nutrients are homogeneously set to their plantless equilibrium level H0 = p/e and N0 = Nin/r, respectively[19]. The system required about 150 yr to reach a quasi-steady-state pattern, as displayed in Supplementary Movie 1.

The pattern analysis was performed on vascular plant and Sphagnum biomass data, to infer “ridge” and “hollow” components of the landscape. The patterns were binarized, with a threshold at 500 and 200 gm−2 for vascular plants and Sphagnum, respectively. The patterns, the nutrient inputs ensuring the transitions between pattern types, and the critical exponents close to the percolation transition are quite insensitive to these thresholds, over a wide range of values (100–1400 g m−2 and 50–600 g m−2, respectively for vascular plants and Sphagnum). In contrast, the absolute values of the filling factors varied with these thresholds. We then detected each cluster of the considered species as a connected component using the Matlab function bwconncomp and determined their characteristics (size, surface, shape, etc.) using the Matlab function regionprops. The number of clusters of both vascular plants and Sphagnum, displayed as a function of nutrient input Nin, was fitted with Gaussian functions. As detailed in the Results section above, the peak of each Gaussian function defines the transition from approximately circular to elongated clusters.

We defined percolating clusters as clusters continuously connecting two opposite sides of the simulation grid. Due to the periodic boundary conditions, this definition does not exactly match the standard one: percolating clusters cover a full turn around either the small or the large diameter of the torus defined by the periodic boundary conditions. The jump from percolating Sphagnum to percolating vascular plants defines a threshold in Nin. As for other experimental systems like filamentation patterns in high-power laser beams[30] or liquid crystals[38], the approach of the percolation point was characterized with a critical exponent γ, determined by fitting the average cluster size (excluding the percolating cluster) against the filling factor (i.e., the fraction of the surface covered by the considered species) with a power law[38].