Quantitative analysis of self-organized patterns in ombrotrophic peatlands

TL;DR

This paper uses numerical modeling to analyze and classify self-organized vegetation patterns in peatlands, focusing on how nutrient levels influence pattern transitions and stability under climate stress.

Contribution

It introduces a systematic quantitative method to classify peatland vegetation patterns based on statistical analysis of clusters, linking pattern transitions to nutrient availability.

Findings

Identified transition from Sphagnum to vascular plant patterns with increased nutrients.

Developed a quantitative framework for pattern stability analysis under drought conditions.

Characterized pattern stability and transitions in peatlands under climate stress.

Abstract

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 \textit{Sphagnum} mosses and vascular plants. In particular, we define the transition from \textit{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.

Figures

Tables

Table 1. Table 1 : Parameter and variables symbols, units and default values  [ 10 , 20 ] . Superscript ∗ denotes values affected by climate. g B subscript g B \mathrm{g}_{\mathrm{B}} , g S subscript g S \mathrm{g}_{\mathrm{S}} , and g N subscript g N \mathrm{g}_{\mathrm{N}} respectively denote grams of vascular biomass, Sphagnum biomass, and nutrients.

Symbol	Interpretation	Unit	Default value
Variables
$B$	Vascular plant piomass	$g_B{m}^{-2}$	Variable
$S$	Sphagnum biomass	$g_S{m}^{-2}$	Variable
$H$	Hydraulic head	$m$	Variable
$N$	Nutrient availability	$g_N{m}^{-2}$	Variable
$N_{\text{c}}$	Nutrient concentration	$g_N{m}^{-3}$	Variable
water-stress function: $f\left(H\right)=\frac{H-z-h_2}{h_1-h_2}$ for $h_1\ge H-z\ge h_2$, $f\left(H\le h_2+z\right)=0$, $f\left(H\ge h_1+z\right)=1$
$h$	Pressure head	$m$	Variable
$h_1$	Pressure head below which soil water stress occurs	$m$	0
$h_2$	Root depth	$m$	-0.5
$z$	Elevation head	$m$	1
Vascular biomass equation
$r_{\text{B}}$	Vascular plant growth parameter	$m^3{g}_N^{-1}{\mathrm{yr}}^{-1}$	0.2
$d$	Vascular plant recycling parameter	${\mathrm{yr}}^{-1}$	0.1
$b$	Vascular plant loss parameter	${\mathrm{yr}}^{-1}$	0.2
$D_{\text{B}}$	Vascular plant diffusion parameter	$m^2{\mathrm{yr}}^{-1}$	2
$c_{\text{BS}}$	Sphagnum inhibiting coefficient on vascular biomass	$g_B{g}_S^{-1}{\mathrm{yr}}^{-1}$	0.1
$k_{\text{BS}}$	Feedback coefficient for vascular plants	$g_B{g}_S^{-1}{\mathrm{yr}}^{-1}$	0.5
$H_{\text{B}}$	Half-saturation density for vascular plants	$g_B{m}^{-2}$	800
Sphagnum biomass equation
$r_{\text{S}}$	Sphagnum growth parameter	${\mathrm{yr}}^{-1}$	0.2
$S_{\text{max}}$	Maximum density of Sphagnum mosses	$g_S{m}^{-2}$	800
$k_{\text{SB}}$	Feedback coefficient for Sphagnum	$g_S{g}_B^{-1}{\mathrm{yr}}^{-1}$	0.1
$H_{\text{S}}$	Half-saturation density for vascular plants	$g_S{m}^{-2}$	300
$D_{\text{S}}$	Diffusion coefficient for Sphagnum	$m^2{\mathrm{yr}}^{-1}$	0.2
Hydraulic head equation
$p$	Precipitation	${\mathrm{myr}}^{-1}$	0.5∗
$t_{\text{v}}$	Transpiration parameter	$m^3{g}_B^{-1}{\mathrm{yr}}^{-1}$	0.005∗
$e$	Evaporation parameter	${\mathrm{myr}}^{-1}$	0.3∗
$k$	Hydraulic conductivity	${\mathrm{myr}}^{-1}$	500
$\theta$	Soil porosity	Dimensionless	0.7
Nutrient availability equation
$N_{\text{in}}$	Nutrient input	$g_N{m}^{-2}{\mathrm{yr}}^{-1}$	0–5
$u$	Vascular plant uptake parameter	$m^3{g}_B^{-1}{\mathrm{yr}}^{-1}$	0.002
$r$	Nutrient loss parameter	${\mathrm{yr}}^{-1}$	0.1
$D_{\text{N}}$	Diffusion coefficient for nutrients	$m^2{\mathrm{yr}}^{-1}$	10

Equations

$$N=N_{\text{c}}H\theta$$

$$\frac{\partial B}{\partial t}=r_{\text{B}}{N}_{\text{c}}Bf(H)-(d+b)B-c_{\text{BS}}S-\frac{k_{\text{BS}}SB}{H_B+B}+D_{\text{B}}\left(\frac{{\partial }^{2}B}{\partial x^2}+\frac{{\partial }^{2}B}{\partial y^2}\right)$$

$$\frac{\partial N}{\partial t}=N_{\textrm{in}}-uN_{\textrm{c}}Bf(H)+\frac{du}{g}B-rN+\theta D_{\textrm{N}}H\left(\frac{\partial^{2}N_{\textrm{c}}}{\partial x^{2}}+\frac{\partial^{2}N_{\textrm{c}}}{\partial y^{2}}\right)+kH\biggl{[}\frac{\partial}{\partial x}\left(N_{\textrm{c}}\frac{\partial H}{\partial x}\right)+\frac{\partial}{\partial y}\left(N_{\textrm{c}}\frac{\partial H}{\partial y}\right)\biggl{]}$$

$$\frac{\partial H}{\partial t}=\frac{p}{\theta}-\frac{e}{\theta}f(H)-\frac{{t_{\textrm{v}}}}{\theta}f(H)B+\frac{k}{\theta}\biggl{[}\frac{\partial}{\partial x}\left(H\frac{\partial H}{\partial x}\right)+\frac{\partial}{\partial y}\left(H\frac{\partial H}{\partial y}\right)\biggl{]}$$

$$\frac{\partial S}{\partial t}=r_{\text{s}}S\left(1-\frac{S}{S_{\text{max}}}\right)-\frac{k_{\text{SB}}SB}{H_{\text{S}}+S}+D_{\text{S}}\left(\frac{{\partial }^{2}S}{\partial x^2}+\frac{{\partial }^{2}S}{\partial y^2}\right)$$

$$\mathcal{A}\propto |p-p_{\text{c}}{|}^{-\gamma }$$

Full text

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 $4\,000\,000\text{\,}\mathrm{k}\mathrm{m}^{2}$ worldwide, they constitute an essential carbon well storing $270550\text{\,}\mathrm{P}\mathrm{g}\mathrm{C}$ [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, 7, 8, 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 [13, 14, 15, 4]. Based on two-dimensional modelling, several processes have been proposed to create or influence these patterns, including peat accumulation, water ponding, and nutrient accumulation [16, 9, 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, 22, 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 (Figure 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.

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 $B(x,y,t)$ and $S(x,y,t)$, hydraulic head $H(x,y,t)$, and nutrient availability $N(x,y,t)$, the latter being linked to the nutrient concentration per unit surface $N_{\textrm{c}}(x,y,t)$ (both in $\mathrm{g}\mathrm{/}\mathrm{m}^{2}$) by

[TABLE]

where $\theta$ describes the soil porosity. The vascular plant biomass density evolves as:

[TABLE]

where the first term is the growth of the plants, $r_{\textrm{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_{\textrm{in}}$ mostly originates from the atmosphere, via precipitation [24, 1]. 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

[TABLE]

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 $\theta$), the evaporation parameter $e$, and the transpiration parameter $t_{\textrm{v}}$:

[TABLE]

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

[TABLE]

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\times\textrm{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\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$\times$3.75\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/} size centred on $90\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ East, $61.23\text{\,}\mathrm{\SIUnitSymbolDegree}\text{/}$ 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_{\textrm{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 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_{\textrm{in}}$: (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 $400\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$. 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).

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_{\textrm{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_{\textrm{in}}$ (Figure 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_{\textrm{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), 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_{\textrm{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 $\mathcal{A}$ (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 (Figure 3b) in the vicinity of the transition. This transition occurs at the critical filling factor $p_{\textrm{c}}$ ($p_{\textrm{c,V}}=0.37$ and $p_{\textrm{c,S}}=0.63$ for binarization thresholds of 500 and $200\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}^{2}$ for vascular and Sphagnum biomass, respectively):

[TABLE]

We obtain $\gamma_{\textrm{V}}=0.83\pm 0.08$ and $\gamma_{\textrm{S}}=0.71\pm 0.15$ for the vascular plant and Sphanum biomass, respectively (Figure 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\approx 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.

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_{\textrm{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 Figure 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 approaches 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 (Figure 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_{\textrm{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_{\textrm{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\times\textrm{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\text{\,}\mathrm{g}\mathrm{m}^{-2}\mathrm{y}\mathrm{r}^{-1}$. The system equilibrates by adapting the vascular plant biomass (Figure 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 (Figure 6c,d). This feedback is faster (at the year-scale) for high nutrient availability ($N_{in}\geq$0.75\text{,}\mathrm{g}\mathrm{m}^{-2}\mathrm{y}\mathrm{r}^{-1}), and slows down when the input nutrient flux is lower (see [Figure 6](#Sx3.F6)). 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 [Figures 2](#Sx3.F2) and [3](#Sx3.F3)). 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 ([Figure 6](#Sx3.F6)a,c) with a certain time-lag. For $N_{in}\geq$0.75\text{\,}\mathrm{g}\mathrm{m}^{-2}\mathrm{y}\mathrm{r}^{-1}, there is enough vascular biomass in the system to respond rapidly (yearly time scale) to variations in precipitation (Figure 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\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$ to $50\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$. 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 $10\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$. Correspondingly, the nutrient intake by vascular plants is reduced, so that nutrients accumulate in the soil during the drought (Figure 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.254.75\text{\,}\mathrm{g}\mathrm{m}^{-2}\mathrm{y}\mathrm{r}^{-1}$, corresponding to the Regimes 2–6 displayed in Figures 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 [34, 22], including non-ombrotrophic peatlands [35], mussel banks [11], coral reefs [36], drylands [37], or savanna [38].

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

Methods

Numerical implementation

Simulations were performed in Matlab on a 128 $\times$ 128 pixel grid with periodic (toroidal) boundary conditions, covering 256 $\times$ $256\text{\,}\mathrm{m}$, using a Newton (Euler) numerical scheme with a temporal step d$t$=${10}^{-4}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$. We started from an initial random distribution of spots bearing $200\text{\,}\mathrm{g}\mathrm{m}^{-2}$ biomass of either vascular plants or Sphagnum mosses, each covering $12.5\text{\,}\%$ of the grid elements. We checked that the results are independent from the particular initial random pattern. The hydraulic head and nutrients are homogeneously set to their plantless equilibrium level $H_{0}=p/e$ and $N_{0}=N_{\textrm{in}}/r$, respectively [19]. The system required about $150\text{\,}\mathrm{y}\mathrm{r}$ 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\text{\,}\mathrm{g}\mathrm{m}^{-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\text{\,}\mathrm{g}\mathrm{m}^{-2}$ and 50 – $600\text{\,}\mathrm{g}\mathrm{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_{\textrm{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_{\textrm{in}}$. As for other experimental systems like filamentation patterns in high-power laser beams [30] or liquid crystals [39], the approach of the percolation point was characterized with a critical exponent $\gamma$, 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 [39].

Acknowledgements

This work was supported by the FNS (grants 200021_155970 and 200020_175697). We gratefully acknowledge fruitful discussions with S. Goyette, G. Rohat, and technical support by M. Moret. Part of the computations were performed at University of Geneva on the Baobab cluster.

Author contributions statement

JK designed the study; CB developed the model and performed the simulations; CB, MB and JK analyzed the results; All authors reviewed the manuscript.

Competing Interests

The authors declare no competing interests

Data availability statement

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

Supplementary Materials

Supplementary Movie 1: Temporal evolution of vascular and Sphagnum mosses biomass, nutrient availability and hydraulic head over 400 years in an ombrotrophic peatland. Initial conditions and parameters as specified in the Methods section and in Table 1. Nutrient input N_{in}=$$2.5\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}^{2}\mathrm{/}\mathrm{y}\mathrm{r}. Each panel represents a square 256-meters wide area, with $2\text{\,}\mathrm{m}$ resolution.

Supplementary Movie 2: Temporal evolution of vascular plant patterns in an ombrotrophic peatland area of 256 by 256 m, after the switch to the $2\times$ CO2 scenario of GCMII, for the values of $N_{\textrm{in}}$ corresponding to the different patterns of Figure 3: $N_{\textrm{in}}$ = 0.25, 1.25, 2.25, 3, 3.875, and $4.5\text{\,}\mathrm{g}\mathrm{/}\mathrm{m}^{2}\mathrm{/}\mathrm{y}\mathrm{r}$. Years 415–424 correspond to a drought, with a drop in precipitation to $50\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{y}\mathrm{r}$.