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.

where θ describes the soil porosity. The vascular plant biomass density evolves as: where the first term is the growth of the plants, r B 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 h 1 + 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. The growth of vascular plants is regulated by their access to nutrients. In an ombrotrophic peatland, the nutrient external input N in 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 as where 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 t v : 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 by where 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][27][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 × CO 2 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 N in 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 N in :   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][22][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, N in .
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 N in (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 N in . 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.
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),  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 in = 0.25, 1.25, 2.25, 3, 3.875, and 4.5 g/m 2 /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 N in = 2.5 g/m 2 /yr. 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 S max , 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) (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 ≈ . 43/18 2 34 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.  Table 1  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 t v 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][14][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 N in , 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 t v 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.
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 × CO 2 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)  , 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. 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.
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.  Table 1

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 H 0 = p/e and N 0 = N in /r, respectively 19 . The system required about 150 yr to reach a quasi-steady-state pattern, as displayed in Supplementary Movie 1. pattern analysis. 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 N in , 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 N in . 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 .

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.