Modelling the habitat preference of two key Sphagnum species in a poor fen as controlled by capitulum water content

Current peatland models generally treat vegetation as static, although plant community structure is known to alter as a response to environmental change. Because the vegetation structure and ecosystem functioning are tightly linked, realistic projections of peatland response to climate change require the inclusion of vegetation dynamics in ecosystem models. In peatlands, Sphagnum mosses are key engineers. Moss community composition primarily follows habitat moisture conditions. The known species habitat preference along the prevailing moisture gradient might not directly serve as a reliable predictor for future species compositions, as water table fluctuation is likely to increase. Hence, modelling the mechanisms that control the habitat preference of Sphagna is a good first step for modelling community dynamics in peatlands. In this study, we developed the Peatland Moss Simulator (PMS), which simulates the community dynamics of the peatland moss layer. PMS is a process-based model that employs a stochastic, individual-based approach for simulating competition within the peatland moss layer based on species differences in functional traits. At the shoot-level, growth and competition were driven by net photosynthesis, which was regulated by hydrological processes via the capitulum water content. The model was tested by predicting the habitat preferences of Sphagnum magellanicum and Sphagnum fallax – two key species representing dry (hummock) and wet (lawn) habitats in a poor fen peatland (Lakkasuo, Finland). PMS successfully captured the habitat preferences of the two Sphagnum species based on observed variations in trait properties. Our model simulation further showed that the validity of PMS depended on the interspecific differences in the capitulum water content being correctly specified. Neglecting the water content differences led to the failure of PMS to predict the habitat preferences of the species in stochastic simulations. Our work highlights the importance of the capitulum water content with respect to the dynamics and carbon functioning of Sphagnum communities in peatland ecosystems. Thus, studies of peatland responses to changing environmental conditions need to include capitulum water processes as a control on moss community dynamics. Our PMS model could be used as an elemental design for the future development of dynamic vegetation models for peatland ecosystems.

1 Introduction

Peatlands have important roles in the global carbon cycle, as they store about 30 % of the world's soil carbon (Gorham, 1991; Hugelius et al., 2013). Environmental changes, such as climate warming and land use changes, are expected to impact the carbon functioning of peatland ecosystems (Tahvanainen, 2011). Predicting the functioning of peatlands under environmental change conditions requires models to quantify the interactions among ecohydrological, ecophysiological, and biogeochemical processes. These processes are known to be strongly regulated by vegetation (Riutta et al., 2007; Wu and Roulet, 2014), which can change over decadal timescales under changing hydrological conditions (Tahvanainen, 2011). Peatland models have generally considered vegetation structure in an unrealistic manner: as a static component (e.g. Frolking et al., 2002; Wania et al., 2009). The recent regional-scale peatland model developed by Chaudhary et al. (2017) includes dynamic vegetation shifts among a single moss plant functional type (PFT) and four vascular PFTs; however, in order to support realistic predictions of peatland functioning and global biogeochemical cycles, the mechanisms that drive changes in moss community structure need to be identified and integrated with ecosystem processes.

A major fraction of peatland biomass is formed by Sphagnum mosses (Hayward and Clymo, 1983; Vitt, 2000). Although individual Sphagnum species often have narrow habitat niches (Johnson et al., 2015), different Sphagnum species replace each other along the water table gradient; therefore, as a genus, Sphagnum species are spread across a wide range of water table conditions (Rydin and McDonald, 1985; Andrus, 1986; Rydin, 1993; Laine et al., 2009). The species composition of the Sphagnum community strongly affects ecosystem processes such as carbon sequestration and peat formation through interspecific variability in species traits including the photosynthetic potential and litter quality (Clymo, 1970; O'Neill, 2000; Vitt, 2000; Turetsky, 2003). The Sphagnum biomass and litter production gradually raises the moss carpet, which feeds back into the species composition (Robroek et al., 2009). Hence, modelling the moss community dynamics is fundamental for predicting temporal changes in peatland vegetation. As the distribution of Sphagnum species primarily follows the variability in the peatland water table (Andrus, 1986; Väliranta et al., 2007), modelling the habitat preference of Sphagnum species along a moisture gradient could be a good first step in predicting moss community dynamics (Blois et al., 2013).

For a given Sphagnum species, the optimal habitat represents the environmental conditions for it to achieve higher rates of net photosynthesis and shoot elongation than its peers (Titus and Wagner, 1984; Rydin and McDonald, 1985; Rydin, 1997; Robroek et al., 2007a; Keuper et al., 2011). The capitulum water content and water storage, which are determined by the balance between the evaporative loss and water gains from capillary rise and precipitation, represent one of the most important controls on net photosynthesis (Titus and Wagner, 1984; Murray et al., 1989; Van Gaalen et al., 2007; Robroek et al., 2009). To quantify the water processes in mosses, hydrological models have been developed to simulate the water movement between the moss carpet and the peat underneath, as regulated by the variations in meteorological conditions and energy balance (Price, 2008; Price and Whittington, 2010). On the other hand, experimental work has addressed the species-specific responses of net photosynthesis to changes in the capitulum water content (Titus and Wagner, 1984; Hájek and Beckett, 2008; Schipperges and Rydin, 1998) and light intensity (Rice et al., 2008; Laine et al., 2011; Bengtsson et al., 2016). Net photosynthesis and hydrological processes are linked via capitulum water retention, which controls the response of the capitulum water content to water potential changes (Jassey and Signarbieux, 2019). However, these mechanisms have not been integrated with ecosystem processes in model simulations.

Along with the capitulum water processes, modelling the habitat preference of Sphagna requires the quantification of the competition among mosses, i.e. the “race for space” (Rydin, 1993, 1997; Robroek et al., 2007a; Keuper et al., 2011): Sphagnum shoots can form new capitula and spread laterally if there is space available. This reduces or eliminates the light source for any plant that is buried by its peers (Robroek et al., 2009). As the competition occurs between neighbouring shoots, its modelling requires downscaling interlinked water–energy processes from the ecosystem to the shoot level. Thus, Sphagnum competition needs to be modelled as spatial processes, considering the fact that spatial coexistence and the variations in functional traits among shoot individuals may impact the community dynamics (Bolker et al., 2003; Amarasekare, 2003). However, coexistence generally relies on simple coefficients to describe the interactions among individuals (e.g. Czárán and Iwasa, 1998; Anderson and Neuhauser, 2000; Gassmann et al., 2003; Boulangeat et al., 2018) and is consequently decoupled from environmental fluctuation or the stochasticity of biophysiological processes.

This study aims to develop and test a model, the Peatland Moss Simulator (PMS), to simulate community dynamics within the peatland moss layer that results in realistic habitat preference of Sphagnum species along a moisture gradient. In PMS, community dynamics is driven by Sphagnum photosynthesis. Photosynthesis, in turn, is regulated by capitulum water retention through the capitulum moisture content. Therefore, we hypothesize that the water retention of the capitula is the mechanism driving moss community dynamics. We test the model validity using data from an experiment based on two Sphagnum species that have different positions along the moisture gradient in the same peatland site. If our hypothesis holds, the model will (1) correctly predict the competitiveness of the two species in wet and dry habitats and (2) fail to predict competitiveness if the capitulum water retention and water content of the two species are not correctly specified.

2 Materials and methods

2.1 Study site

The peatland site being modelled is located in Lakkasuo, Orivesi, Finland (61^∘47^′ N; 24^∘18^′ E). The site is a poor fen fed by mineral inflows from a nearby esker (Laine et al., 2004). Most of the site is formed by lawns dominated by Sphagnum recurvum complex (Sphagnum fallax accompanied by Sphagnum flexuosum and Sphagnum angustifolium) and Sphagnum papillosum. Less than 10 % of the surface is occupied by hummocks, inhabited by Sphagnum magellanicum and Sphagnum fuscum, which are 15–25 cm higher than the lawn surfaces. Both microforms are covered by continuous Sphagnum carpet with a sparse vascular plant cover (12 % Carex cover on average), which spread homogeneously over the topography. The annual mean water table was 15.6±5.0 cm deep at the lawn surface (Kokkonen et al., 2019). More information about the site can be found in Kokkonen et al. (2019).

2.2 Model outline

The Peatland Moss Simulator (PMS) is a process-based, stochastic model that simulates the temporal dynamics of a Sphagnum community as driven by variations in precipitation, irradiation, and energy flow with individual-based interactions (Fig. 1). In PMS, the studied ecosystem is seen as a dual-column system consisting of hydrologically connected habitats of hummocks and lawns (community environment in Fig. 1). For each habitat type, the community area is downscaled to two-dimensional cells representing the scale of individual shoots (i.e. 1 cm²). Each grid cell can be occupied by one capitulum from a single Sphagnum species. The community dynamics, i.e. the changes in species abundances, are driven by the growth and competition of Sphagnum shoots at the grid-cell level (Module I in Fig. 1). These processes are regulated by the grid-cell-specific conditions of water and energy (Module II in Fig. 1), which are derived from the community environment (Module III in Fig. 1).

Figure 1 Framework of the Peatland Moss Simulator (PMS).

In this study, we focused on developing modules I and II (Sect. 2.3) and employed an available soil–vegetation–atmosphere transport (SVAT) model (Gong et al., 2013, 2016) to describe the water–energy processes for Module III (Appendix A). We assumed that the temporal variation in the water table was similar in lawns and hummocks and that the hummock–lawn differences in the water table (dWT in Fig. 1) followed their difference in surface elevations (Wilson, 2012). At the grid-cell level, the photosynthesis of capitula drove the biomass growth and elongation of shoots, which led to the competition between adjacent grid cells. The net photosynthesis rate was controlled by the capitulum water content, W_cap, which was defined by the capitulum water retention in relation to the water potential, h (Sect. 2.4). The values for functional traits that regulate the growth and competition processes were randomly selected within their normal distribution measured in the field (Sect. 2.4). Unknown parameters that related the lateral water flows of the site are estimated using a machine-learning approach (Sect. 2.5). Finally, a Monte Carlo simulation was used to support the analysis of the habitat preferences of Sphagnum species and hypothesis tests (Sect. 2.6). The list of symbols used is given in Appendix E.

Table 1 Species-specific values of morphological and photosynthetic parameters for S. magellanicum and S. fallax. The parameters include the capitulum density (D_S), the capitulum biomass (B_cap), the specific height of the stem (H_spc), the maximal gross photosynthesis rate at 20 ^∘C (P_m20), the respiration rate at 20 ^∘C (R_s20), the half-saturation point of photosynthesis (α_PPFD), and the polynomial coefficients (a_W0, a_W1, and a_W2) for the responses of net photosynthesis to capitulum water content. Parameter values are given as the mean ± standard deviation.

^∗ The parameter was used in the linear models predicting the log₁₀-transformed capitulum water potential (h) and the bulk resistance (r_bulk) for S. fallax and S. magellanicum. The capitulum density and photosynthetic parameter values measured here are well within the range of those reported in the literature for these species (McCarter and Price, 2014; Laing et al., 2014; Bengtsson et al., 2016; Korrensalo et al., 2016).

2.3 Model development

2.3.1 Calculating shoot growth and competition of Sphagnum mosses (Module I)

Calculation of Sphagnum growth

To model the grid-cell biomass production and height increment, we assumed that capitula were the main parts of the shoots responsible for photosynthesis and the production of new tissues, instead of the stem sections underneath. We employed a hyperbolic light-saturation function (Larcher, 2003) to calculate the net photosynthesis, which was parameterized based on empirical measurements made from the target species collected from the study site (see Appendix B for materials and methods):

$$\begin{array}{}\text{(1)}& A_{\mathrm{20}}=\left({\displaystyle \frac{P_{\mathrm{m}\mathrm{20}}\times \text{PPFD}}{{\mathit{\alpha }}_{\text{PPFD}}+\text{PPFD}}}-R_{\mathrm{s}\mathrm{20}}\right)\times B_{\text{cap}},\end{array}$$

where subscript 20 denotes the variable value measured at 20 ^∘C; R_s is the mass-based respiration rate (µmol g⁻¹ s⁻¹); P_m is the mass-based rate of maximal gross photosynthesis (µmol g⁻¹ s⁻¹); PPFD is the photosynthetic photon flux density (µmol m⁻² s⁻¹); B_cap is the capitulum biomass; and α_PPFD is the half-saturation point (µmol m⁻² s⁻¹) for photosynthesis.

By adding multipliers for the capitula water content (f_W) and temperature (f_T) to Eq. (1), the net photosynthesis rate A (µmol m⁻² s⁻¹) was calculated as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& A=\left[{\displaystyle \frac{{P_{\mathrm{m}}}_{\mathrm{20}}\times \text{PPFD}}{{\mathit{\alpha }}_{\text{PPFD}}+\text{PPFD}}}{f}_{\text{T}}\left(T\right)-R_{\mathrm{s}\mathrm{20}}{f}_{\mathrm{R}}\left(T\right)\right]\times B_{\text{cap}}\\ & \times f_{\text{W}}\left(W_{\text{cap}}\right),\end{array}\end{array}$$

where f_W(W_cap) describes the responses of A to the capitulum water content, W_cap; and f_T(T) describes the responses of P_m to the capitulum temperature T (Korrensalo et al., 2017). f_W(W_cap) was estimated based on the empirical measurements (see Sect. 2.4 and Appendix B). The temperature response f_R(T) is a Q₁₀ function that describes the temperature sensitivity of R_s (Frolking et al., 2002):

$$\begin{array}{}\text{(3)}& f_{\mathrm{R}}\left(T\right)=Q_{\mathrm{10}}^{\left(T-T_{\mathrm{opt}}\right)/\mathrm{10}},\end{array}$$

where Q₁₀ is the sensitivity coefficient; T is the capitulum temperature (^∘C); and T_opt (20 ^∘C) is the reference temperature of respiration.

The response of A to W_cap (f_W(W_cap); Eq. 2) was described as a second-order polynomial function:

$$\begin{array}{}\text{(4)}& f_{\text{W}}\left(W_{\text{cap}}\right)=a_{\text{W0}}+a_{\text{W1}}\times W_{\text{cap}}+a_{\text{W2}}\times W_{\text{cap}}^{\mathrm{2}},\end{array}$$

where a_W0, a_W1, and a_W2 are coefficients.

Plants can store carbohydrates as nonstructural carbon (NSC, e.g. starch and soluble sugar) to support fast growth in spring or during post-stress periods, such as after drought events (Smirnoff, 1992; Martínez-Vilalta et al., 2016; Hartmann and Trumbore, 2016). We linked the production of shoot biomass to the immobilization of NSC storage (modified from Eq. 10 in Asaeda and Karunaratne, 2000). The change in NSC storage depends on the balance between net photosynthesis and immobilization:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& M_{\mathrm{B}}=s_{\text{imm}}\times \mathrm{NSC}\times k_{\text{imm}}{\mathit{\alpha }}_{\text{imm}}^{T-\mathrm{20}}\\ & \partial \text{NSC}/\partial t=A-M_{\mathrm{B}},\text{NSC}\mathit{ϵ}\left[\mathrm{0},{\text{NSC}}_{\text{max}}\right],\end{array}\end{array}$$

where M_B is the immobilized NSC for biomass production during a time step (g); k_imm is the specific immobilization rate (g g⁻¹) (Asaeda and Karunaratne, 2000); α_imm is the temperature constant; s_imm is the multiplier for the temperature threshold, where s_imm=1 when T>5 ^∘C but s_imm=0 if T≤5 ^∘C; and NSC_max is the maximal NSC concentration in Sphagnum biomass (Turetsky et al., 2008). The timing of growth is controlled by a temperature threshold and NSC availability. Growth occurs when T>5 ^∘C and NSC is above zero. The dynamics of NSC storage is related to the water content (WC) through net photosynthesis.

The increase in shoot biomass drove the shoot elongation:

$$\begin{array}{}\text{(6)}& \partial H_{\mathrm{c}}/\partial t={\displaystyle \frac{M_{\mathrm{B}}}{H_{\text{spc}}{S}_{\text{c}}}},\end{array}$$

where H_c is the shoot height (cm); H_spc is the biomass density of Sphagnum stems (g m⁻² cm⁻¹); and S_c is the area of a cell (m²).

Calculation of Sphagnum competition and community dynamics

To simulate the competition among Sphagnum shoots, we first compared H_c of each grid cell (source grid cell, i.e. grid cell a in Fig. 1) to its four neighbouring cells and marked the one with the lowest position (e.g. grid cell b in Fig. 1) as the target of spreading. The spreading of shoots from a source to a target grid cell occurred when the following criteria were fulfilled: (i) the height difference between the source and target grid cells exceeded a threshold value; (ii) NSC accumulation in the source grid cell was large enough to support the growth of new capitula in the target grid cell; and (iii) the capitula in the source grid cell could split once per year at most.

The height difference threshold in rule (i) was set equal to the mean diameter of the capitula in the source cell, based on the assumption that the shape of a capitulum was spherical. When shoots spread, the species type and model parameters in the target grid cell were overwritten by those in the source grid cell, assuming the mortality of shoots originally in the target cell. During the spreading, the NSC storage was transferred from the source cell to the target cell to form new capitula. In cases where spreading did not take place, the establishment of new shoots from spores could maintain the continuity of the Sphagnum carpet at the site. During the establishment from spores, which was rare and occurred during the first years of simulation, the traits of Sphagnum species were randomized within their normal distribution measured in the field.

2.3.2 Calculating the grid-cell-level dynamics of environmental factors (Module II)

Module II computes grid-cell values of W_cap, PPFD, and T for Module I. The cell-level PPFD and T were assumed to be equal to the community means, which were solved using the SVAT scheme in Module III (Appendix A). The community-level evaporation rate (E) was partitioned to the cell level (E_i) as follows:

$$\begin{array}{}\text{(7)}& E_i=E\times \left({\displaystyle \frac{{\text{Sv}}_i}{r_{\text{bulk},i}}}\right)/\sum \left({\displaystyle \frac{{\text{Sv}}_i}{r_{\text{bulk},i}}}\right),\end{array}$$

where r_bulk,i is the bulk surface resistance of cell i, which is as a function ($r_{\text{bulk},i}=f_{\mathrm{r}}\left(h_i\right)$) of the grid-cell-based water potential h_i, capitulum biomass (B_cap), and shoot density (D_S) based on the empirical measurements (Appendix B). Sv_i was the evaporative area, which was related to the height differences among adjacent grid cells:

$$\begin{array}{}\text{(8)}& {\text{Sv}}_i={\text{Sc}}_i+l_{\mathrm{c}}{\sum }_j\left(H_{{\mathrm{c}}_i}-H_{{\mathrm{c}}_j}\right),\end{array}$$

where l_c is the width of a grid cell (cm); and subscript j denotes the four nearest neighbouring grid cells. In this way, changes in the height difference between the neighbouring shoots feed back to affect the water conditions of the grid cells via an alteration of the evaporative surface area.

The grid-cell-level changes in the capitula water potential (h_i) were driven by the balance between the evaporation (E_i) and the upward capillary flow to capitula:

$$\begin{array}{}\text{(9)}& \partial h_i={\displaystyle \frac{K_{\mathrm{m}}}{C_i}}\left[{\displaystyle \frac{\left(h_i-h_{\mathrm{m}}\right)}{\mathrm{0.5}{z}_{\mathrm{m}}}}-\mathrm{1}-E_i\right],\end{array}$$

where h_m is the water potential of the living moss layer, solved in Module III (Appendix A); z_m is the thickness of the living moss layer (z_m=5 cm); K_m is the hydraulic conductivity of the moss layer and is set to be the same for each grid cell; and C_i is the cell-level specific water uptake capacity ($C_i=\partial W_{\text{cap},i}/\partial h_i$). $\partial W_{\text{cap},i}/\partial h_i$ could be derived from the capitulum water retention function h_i=f_h(W_cap). W_cap, which affects the calculation of net photosynthesis through f_W(W_cap), can then be estimated from h_i using the f_h(W_cap) function (Table B2, Table B3).

2.4 Model parameterization

2.4.1 Selection of Sphagnum species

We chose S. fallax and S. magellanicum, which form 63 % of total plant cover at the study site at Lakkasuo (Kokkonen et al., 2019), as the target species representing the lawn and hummock habitats respectively. These species share a similar niche along soil pH and nutrient richness gradients (Wojtuń et al., 2003), but they are discriminated by their water table level preferences (Laine et al., 2004): while S. fallax is commonly found close to the water table (Wojtuń et al., 2003), S. magellanicum can occur along a wider range of the wet–dry gradient, from intermediately wet lawns to dry hummocks (Rice et al., 2008; Kyrkjeeide, et al., 2016; Korresalo et al., 2017). Thus, the transition from S. fallax to S. magellanicum along the wet–dry gradient indicates the decreasing competitiveness of S. fallax with S. magellanicum with a lowering water table.

2.4.2 Parameterization of morphological traits, net photosynthesis, and capitulum water retention

We empirically quantified the morphological traits of the capitulum density (D_S; shoots cm⁻²), the biomass of capitula (B_cap; g m⁻²), the biomass density of living stems (H_spc; g cm⁻¹ m⁻²), the net photosynthesis parameters (P_m20, R_s20, and α_PPFD), and the water retention properties, i.e. f_h(W_cap) and f_r(h) (Eqs. 8 and 10), for the two Sphagnum species (see Appendix B for methods). The values (mean ± SD) of the morphological parameters, the photosynthetic parameters, and the polynomial coefficients (a_W0, a_W1, and a_W2; Eq. 3) are listed in Table B1 in Appendix B. For each parameter, a random value was initialized for each cell based on the measured means and SD, assuming that the variation in the parameter values is normally distributed.

Table 2 Parameters values for the SVAT simulations (Module III). The parameters include the saturated hydraulic conductivity (K_sat), the water retention parameters of water retention curves (α and n), the saturated water content (θ_s), the permanent wilting point water content (θ_r), the snow layer surface albedos (a_s, a_l), the thermal conductivity (K_T), the specific heat (C_T), and the maximal nonstructural carbon (NSC) concentration (NSC_max).

^a  The value was calculated from bulk density (ρ_bulk) as ${\mathit{\theta }}_{\mathrm{s}}=\mathrm{97.95}-\mathrm{79.72}{\mathit{\rho }}_{\text{bulk}}$ following Päivänen (1973); ^b The value was calculated as ${\mathit{\theta }}_{\mathrm{r}}=\mathrm{4.3}+\mathrm{67}{\mathit{\rho }}_{\text{bulk}}$ following Weiss et al. (1998).

We noticed that the fitted f_W(W_cap) was meaningful when W_cap was below the optimal water content for photosynthesis $(W_{\text{opt}}=-\mathrm{0.5}{a}_{\text{W1}}/a_{\text{W2}})$. If W_cap>W_opt, photosynthesis decreased linearly with increasing W_cap, as it is limited by the diffusion of CO₂ (Schipperges and Rydin, 1998). In that case, f_W(W_cap) was calculated following Frolking et al. (2002):

$$\begin{array}{}\text{(10)}& f_{\text{W}}\left(W_{\text{cap}}\right)=\mathrm{1}-\mathrm{0.5}{\displaystyle \frac{W_{\text{cap}}-W_{\text{opt}}}{W_{\text{max}}-W_{\text{opt}}}},\end{array}$$

where W_max is the maximum water content of the capitula.

It is known that W_max is around 25–30 g g⁻¹ (e.g. Schipperges and Rydin, 1998) or about 0.31–0.37 cm³ cm⁻³ in terms of the volumetric water content (assuming a 75 g m⁻² capitula biomass and a 0.6 cm high capitula layer). This range is broadly lower than the saturated water content of the moss carpet (>0.9 cm³ cm⁻³, McCarter and Price, 2014). Consequently, we used the following equation to convert the volumetric water content to the capitulum water content (RWC), when h_i was higher than the boundary value of −10⁴ cm:

$$\begin{array}{}\text{(11)}& W_{\text{cap}}=\text{min}\left(W_{\text{max}},{\mathit{\theta }}_{\mathrm{m}}/\left(H_{\text{cap}}\times B_{\text{cap}}\times {\mathrm{10}}^{-\mathrm{4}}\right)\right),\end{array}$$

where W_max is the maximum water content and is set to 25 g g⁻¹ for both species; θ_m is the volumetric water content of the moss layer; and H_capis the height of the capitula and is set to 0.6 cm (Hájek and Beckett, 2008).

2.5 Model calibration for lateral water influence

We used a machine-learning approach to estimate the influence of the upstream area on the water balance of the site. The rate of net inflow (I, see Eq. A18 in Appendix A) was described as a function of Julian day (JD), assuming that the inflow was maximum after the spring thaw and then decreased linearly with time:

$$\begin{array}{}\text{(12)}& I_j=\left(a_{\text{N}}\times \text{JD}+b_{\text{N}}\right)\times {K_{\mathrm{s}}}_j\text{JD}>{\text{JD}}_{\text{thaw}},\end{array}$$

where subscript j denotes the peat layers under the water table; K_s is the saturated hydraulic conductivity; JD_thaw is the Julian day that thawing completed; and a_N and b_N are parameters.

We simulated water table changes using climatic scenarios from the Weather Generator (Appendix A). During calibration, the community compositions were set to remain constant, such that S. magellanicum fully occupied the hummock habitat whereas S. fallax fully occupied the lawn habitat. The simulated multiyear means of weekly water table values were compared to the weekly mean water table values observed at the site during the years 2001, 2002, 2004, and 2016. The cost function for the learning process was based on the sum of squared error (SSE) of the simulated water table:

$$\begin{array}{}\text{(13)}& \text{SSE}=\mathrm{\Sigma }{\left({\text{WT}}_{{\mathrm{s}}_k}-{\text{WT}}_{{\mathrm{m}}_k}\right)}^{\mathrm{2}},\end{array}$$

where WT_m is the measured multiyear weekly mean of the water table; WT_s is the simulated multiyear weekly mean of the water table; and subscript k denotes the week of the year when the water table was sampled.

The values of a_N and b_N were estimated using the gradient descent approach (Ruder, 2016) by minimizing the SSE in a_N and b_N in Eq. (12):

$$\begin{array}{}\text{(14)}& X_{\text{N}}\left(j\right):=X_{\text{N}}\left(j\right)-\mathrm{\Gamma }{\displaystyle \frac{\partial \text{SSE}}{\partial X_{\text{N}}\left(j\right)}},\end{array}$$

where Γ is the learning rate (Γ=0.1). Appendix D shows the simulated water table with the calibrated inflow term I compared with the measured values from the site.

2.6 Model-based analysis

First, we examined the ability of the model to capture the preference of S. magellanicum for the hummock environment and the preference of S. fallax for the lawn environment (Test 1). For both species, the probability of occupation was initialized as 50 % in a cell, and the distribution of the species in the communities were randomly patterned. Monte Carlo simulations (40 replicates) were carried out with a time step of 30 min. A simulation length of 15 years was selected based on preliminary studies in order to cover the major interval of change and to ease computational demand. Biomass growth, stem elongation, and the spreading of shoots were simulated on a daily basis. The establishment of new shoots in deactivated cells was calculated at the end of each simulation year. We then assessed if the model could capture the dominance of S. magellanicum in the hummock communities and the dominance of S. fallax in lawn communities. The simulated annual height increments of mosses were compared to the values measured for each community type. To measure moss vertical growth in the field, we deployed 20 cranked wires on S. magellanicum-dominated hummocks and 15 cranked wires on S. fallax-dominated lawns in 2016. Each cranked wire was a piece of metal wire attached to plastic brushes at the side that was anchored into the moss carpet (e.g. Clymo, 1970; Holmgren et al., 2015). Annual vertical growth (dH) was determined by measuring the change in the exposed wire length above the moss surface from the beginning to the end of growing season.

Table 3 The results from Test 2, which addressed the robustness of the model to the uncertainties in a set of parameters. Each parameter was increased or decreased by 10 %. The model was run for S. magellanicum and S. fallax in their preferential habitats. The difference in mean cover between simulations as a function of changed or unchanged parameter values is given; the standard deviations (SD) of the means are given in parentheses. The parameters include the specific immobilization rate (kimm), the maximal nonstructural carbon (NSC) concentration (NSC_max), the hydraulic conductivity of the moss layer (K_m), the hydraulic conductivity of the peat layer (K_h), the water retention parameters of the water retention curves (α and n), the snow layer surface albedo (a_s), and the aerodynamic resistance (r_aero).

Second, we tested the robustness of the model to the uncertainties in a set of parameters (tests 2–4). In Test 2, we focused on parameters that were closely linked to hydrology and growth calculations but were roughly parameterized (e.g. k_imm, r_aero) or adopted as a priori from other studies (e.g. K_sat, α, n, NSC_max; see Table 2). One at a time, each parameter value was adjusted by +10 % or −10. A total of 40 Monte Carlo simulations were run using the same runtime settings as in Test 1. The simulated means of cover were then compared to those calculated without the parameter adjustment.

Tests 3 and 4 were then carried out to test whether the model could correctly predict the competitiveness of the species in dry and wet habitats if the species-specific trends in the capitulum water content were not correctly specified. For both species, we set the values of parameters controlling the water retention (i.e. B_cap and D_S; Appendix B) and the water-stress effects on net photosynthesis (i.e. W_cap; Eq. 4) to be the same as those for S. magellanicum (Test 3) or same as those for S. fallax (Test 4). Our hypothesis would be supported if removing the interspecific differences in RWC responses led to the failure to predict the habitat preferences of the species.

We implemented tests 5 and 6 to investigate the importance of parameters that directly control the species ability to overgrow another species with a more rapid height increment (i.e. P_m20, R_s20, α_PPFD, and H_spec) under lawn and hummock conditions. We eliminated the species differences in the parameter, so that the values were the same as those in S. magellanicum (Test 5) or the same as those in S. fallax (Test 6). The effects of the manipulation were compared with those from tests 3 and 4. For tests 3–6, 80 respective Monte Carlo simulations were run using the set-ups described in Test 1.

Tests 7 and 8 were implemented to separate the effects of photosynthetic water-response parameters from the effects of the capitula water retention. We set the photosynthetic water-response parameters to be the same as those in S. magellanicum (Test 7) or the same as those in S. fallax (Test 8). As our model aimed to couple the environmental fluctuations and stochasticity of ecosystem processes, we further tested the model responses to the absence of environmental fluctuations (Test 9) or the absence of stochasticity in model parameters (Test 10). In Test 9, monthly mean values of meteorological variables were used to drive the model simulation. In Test 10, we removed the stochasticity of model parameters and assigned an average value to each parameter of grid cells. For tests 7–10, 40 respective Monte Carlo simulations were run using the set-ups described in Test 1.

3 Results

3.1 Simulating the habitat preferences of Sphagnum species as affected by the capitulum water content traits

Test 1 demonstrated the ability of the model to capture the preference of S. magellanicum for the hummock environment and the preference of S. fallax for the lawn environment (Fig. 2a). The simulated annual changes in species cover were greater in lawn habitats than in hummock habitats during the first 5 simulation years. The changes in lawn habitats slowed down around year 10, and the cover of S. fallax plateaued at around 95±2.8 % (mean ± standard error). In contrast, the cover of S. magellanicum on hummocks continued to grow until the end of the simulation and reached 83±3.1 %. In the lawn habitats, the cover of S. fallax increased in all Monte Carlo simulations, and the species occupied all grid cells in 70 % of the simulations. In the hummock habitats, the cover of S. magellanicum increased in 91 % of Monte Carlo simulations, and it formed a monocultural community in 16 % of simulations (Fig. 2b). The vertical growth of Sphagnum mosses was significantly greater in lawn habitats than in hummock habitats (P<0.01). The ranges of simulated vertical growth agreed well with the observed values from field measurement for both species (Fig. 2c).

Figure 2 Testing the ability of PMS to predict the habitat preference of Sphagnum magellanicum and S. fallax (Test 1). The hummock and lawn habitats were differentiated by the water table depth, the surface energy balances, and the capitulum water potential in the model simulations. At the beginning of the simulation, the cover of the two species was set to be equal and it was allowed to develop with time. (a) The annual development of the relative cover (mean and standard error) of the two species in hummock and lawn habitats, (b) the cumulative probability distribution of the cover of the two species at the end of the 15-year period based on 40 Monte Carlo simulations, and (c) the simulated and measured means of the annual vertical growth of Sphagnum surfaces in their natural habitats (hummock and lawn).

Table 4 Result from tests 7–10, which addressed the importance of meteorological fluctuations, the stochasticity of model parameters, and the photosynthetic water response. In Test 7, monthly mean values of meteorological variables were used to drive the model simulation. In Test 8, the stochasticity of model parameters was removed, and the average values were used for parameters at the grid-cell level. In tests 9 and 10, the photosynthetic water-response parameters (i.e. a_W0, a_W1, and a_W2; see Table 1) were set to be the same as those for S. magellanicum (Test 9) or the same as those for S. fallax (Test 10). The mean cover of S. magellanicum on hummocks and the cover of S. fallax on lawns after the 15-year simulation periods are listed in Table 4.

3.2 Testing model robustness

Test 2 addressed the model robustness to the uncertainties in several parameters that were closely linked to hydrology and growth calculations. Modifying most of the parameter values by +10 % or −10 % yielded marginal changes in the mean cover of species in either hummock or hollow habitats (Table 3). Reducing the moss carpet and the peat hydraulic parameter n had stronger impacts on S. fallax cover in hummock habitats than in lawn habitats. Nevertheless, changes in simulated cover that were caused by parameter manipulations were generally smaller than the standard deviations of the means, i.e. fitting into the random variation.

3.3 Testing the controlling role of the capitulum water content on community dynamics

In tests 3 and 4, the model incorrectly predicted the competitiveness of the two test species when the interspecific differences in the capitulum water content were eliminated. In both tests, S. fallax became dominant in all habitats. The use of water-response characteristics from S. magellanicum for both species (Test 3) led to the faster development of S. fallax cover and higher coverage at the end of simulation (Fig. 3a) compared with the simulation results where the water-response characteristics from S. fallax were used for both species (Test 4, Fig. 3b). The pattern was more pronounced in hummock than in lawn habitats.

Figure 3 Testing the importance of capitulum water content for the habitat preference of S. magellanicum and S. fallax. The development of the relative cover (mean and standard error) were simulated in hummock and lawn habitats over a 15-year time frame for the two species. For both species, parameter values for the capitulum water content, the capitulum biomass (B_cap), and the density (D_S) were set to be the same as those from (a) S. magellanicum (Test 3) or (b) S. fallax (Test 4).

In tests 5 and 6, the species differences in the growth-related parameters were eliminated. However, the model still predicted the dominance of S. fallax and S. magellanicum in lawn and hummock habitats respectively (Fig. 4). The increase in the mean cover of S. magellanicum was especially fast in the hummock habitat in comparison to the results of the unchanged model from Test 1 (Fig. 2a). In lawns, the use of S. fallax growth parameters for both species gave stronger competitiveness to S. magellanicum (Fig. 4b) than using the S. magellanicum parameters (Fig. 4a). In tests 7 and 8, ignoring the interspecific differences in the photosynthetic water-response parameters did not change the simulated habitat preferences of S. fallax and S. magellanicum (Table 4). Using the water-response parameters of S. fallax decreased the mean cover of S. fallax in lawns but increased the cover of S. magellanicum on hummocks. In contrast, using the water-response parameters of S. magellanicum increased the mean cover of S. fallax in lawns but decreased the cover of S. magellanicum on hummocks.

Figure 4 Testing the importance of parameters regulating net photosynthesis and shoot elongation for the habitat preference of S. magellanicum and S. fallax. The annual development of the relative cover (mean and standard error) of the two species was simulated for hummock and lawn habitats over a 15-year time frame. For both species, the parameter values (i.e. P_m20, R_s20, α_PPFD, and H_spec) were set to be the same as those from (a) S. magellanicum (Test 5) or (b) S. fallax (Test 6).

3.4 Testing the effects of environmental fluctuations and the stochasticity of ecosystem processes on community dynamics

In Test 9, the model failed to simulate the preference of S. magellanicum for hummocks (Table 4) if the environmental fluctuation was ignored. However, the simulated cover of S. fallax in lawns was higher compared with the unchanged condition (i.e. Test 1). Using the mean value for each model parameters led to monocultural community, i.e. S. magellanicum occupied 100 % of the hummock area whereas S. fallax took over lawns completely.

4 Discussion

In peatland ecosystems, Sphagna are keystone species distributed primarily along the hydrological gradient (e.g. Andrus, 1986; Rydin, 1986). In a context where substantial change in peatland hydrology is expected under a changing climate in northern areas (e.g. longer snow-free season, lower summer water table, and more frequent droughts), there is a pressing need to understand how peatland plant communities could react and how Sphagnum species could redistribute under habitat changes. In this work, we developed the Peatland Moss Simulator (PMS), a process-based stochastic model, to simulate the competition between S. magellanicum and S. fallax, two key species representing dry (hummock) and wet (lawn) habitats in a poor fen peatland. We empirically showed that these two species differed in characteristics that likely affect their competitiveness along a moisture gradient.

Capitulum water retention for the lawn-preferring species (S. fallax) was weaker than that for the hummock-preferring species (S. magellanicum). Compared with S. magellanicum, the capitula of S. fallax held less water at saturation and the water content decreased more rapidly with dropping water potential. Hence, S. fallax would dry faster than S. magellanicum under the same rate of water loss. Moreover, the water content in S. fallax capitula was less resistant to evaporation. These differences indicated that it is harder for S. fallax capitula to buffer evaporative water loss and, therefore, avoid or delay desiccation. Similar differences between hummock and hollow species have been in previous studies (Titus and Wagner, 1984; Rydin and McDonald, 1985). In addition, the net photosynthesis of S. fallax is more sensitive to changes in the capitulum water content than S. magellanicum, as seen by the steeper decline in photosynthesis with decreasing water content (Fig. B2c). Consequently, the growth of S. fallax is more likely to be slowed down by dry periods, during which the capillary water cannot fully compensate for the evaporative loss (Robroek et al., 2007b), making it less competitive in habitats prone to desiccation.

The PMS successfully captured the habitat preferences of the two Sphagnum species (Test 1): starting from a mixed community with equal probabilities for both species, the lawn habitats with the shallower water table were eventually dominated by the typical lawn species S. fallax, whereas the hummock habitats, which were 15 cm higher than the lawn surface, were taken over by S. magellanicum. The low final cover of S. magellanicum simulated in lawn habitats agreed well with field observations from our study site, where S. magellanicum cover was less than 1 % in lawns (Kokkonen et al., 2019). Conversely, S. fallax was outcompeted by S. magellanicum in the hummock habitats. This result is consistent with previous findings that hollow-preferring Sphagna are less likely to survive in hummock environments with greater drought pressure (see Rydin, 1985; Rydin et al., 2006; Johnson et al., 2015). The simulated annual height increments of mosses also agreed well with the observed values for both habitat types. Our simulation for lawn habitat shows that the looser stem structure of S. fallax allows it to allocate more of its produced biomass into vertical growth and, therefore, overgrow S. magellanicum, in which new biomass forms a compact stem that is packed with thick fascicles. This finding indicates that PMS can capture the key mechanisms controlling the growth and competitive interactions of the Sphagnum species.

Parameter sensitivity testing showed the robustness of PMS regarding the uncertainties in parameterization, as the simulated changes in the mean species cover (under 10 % variation in several parameters) were generally less than the standard deviations of the means. Decreasing the value of the hydraulic parameter n (Table 2, Eq. A5) increased the presence of S. fallax in the hummock habitats. This was expected because n is a scaling factor and its changes are consequently magnified: a lower n value will lead to a higher water content in the unsaturated layers above the water table (van Genuchten, 1980), which allows Sphagna that are adapted to wet environments to survive dry conditions (Hayward and Clymo, 1982; Robroek et al., 2007b; Rice et al., 2008). In contrast, the response of Sphagnum cover to the changes in other hydraulic parameters (i.e. α, n, K_h) was limited in lawn habitats. This could be due to the relatively shallow water table in lawns, which was able to maintain sufficient capillary rise to the moss carpet and capitula. Decreasing the values of the specific immobilization rate (kimm) and maximal NSC concentration in Sphagnum biomass (NSC_max) mainly decreased the cover of S. fallax in lawn habitats, which is consistent with the importance of biomass production to Sphagna in high-moisture environments (e.g. Rice et al., 2008; Laine et al., 2011). In addition, the SVAT model simulations for hummocks and lawns (Module III, Fig. 1) employed the same hydraulic parameter values obtained from S. magellanicum hummocks (McCarter and Price, 2014). For lawns, this could overestimate K_m but underestimate n, as the lawn peat would be less efficient at holding a high water content and generating capillary flow than hummock peat (Robroek et al., 2007b; Branham and Strack, 2014). As the decrease in K_m and increase in n showed counteracting effects on the simulated species covers (Table 3), the biases in the parameterization of K_m and n may not critically impact model performance.

Both our empirical measurements and PMS simulations indicate the importance of the capitulum water content as a mechanism for controlling the moss community dynamics in peatlands. The fact that the Sphagnum niche is defined by two processes has long been hypothesized and experimentally studied. Firstly, dry, high-elevation habitats such as hummocks physically select species with the ability to remain moist (Rydin, 1993). If the interspecific differences in water retention and water-stress effects were correctly specified (Test 1, Fig. 2) our model predicted this phenomena of stronger competitiveness of S. magellanicum with S. fallax in hummock habitats correctly. Alternatively, the model failed to predict the distribution of S. magellanicum on hummocks if these interspecific differences in the water processes were neglected (tests 3 and 4, Fig. 3). During low water table periods in summer, capillary rise may not fully compensate for high evaporation (Robroek et al., 2007b; Nijp et al., 2014). In such circumstances, the capitulum water potential could drop rapidly towards the pressure defined by the relative humidity of air (Hayward and Clymo, 1982). Consequently, the ability of capitula to retain cytoplasmic water is particularly important for the hummock-preferring species, as was also shown by Titus and Wagner (1984).

Secondly, in habitats with a more persistently high moisture content, such as lawns and hollows, interspecific competition becomes important: it is acknowledged that species from such habitats generally have higher growth rates and photosynthetic capacity than hummock species (e.g. Laing et al., 2014; Bengtsson et al., 2016). Our results also agreed with this, as setting the growth-related parameters (i.e. P_m20, R_s20, α_PPFD, and H_spec) of S. magellanicum to be the same as those of S. fallax decreased the S. fallax cover in both hummock and lawn habitats (Test 6, Fig. 4b). However, such changes did not impact the simulated habitat preferences for the tested species. Based on this, the growth-related parameters seem to be less important than water-related parameters. Furthermore, tests 7 and 8 showed that when interspecific differences in the water-stress effects on photosynthesis were removed, the model still predicted correct habitat preferences for S. magellanicum and S. fallax. Therefore, the interspecific differences in the capitulum water retention could be the main determinant of the habitat preferences of the tested species.

There have been growing concerns about the shift of peatland communities from Sphagnum-dominated communities towards higher vascular abundance under a drier and warmer climate (Wullschleger et al., 2014; Munir et al., 2015; Dieleman et al., 2015). Nevertheless, the potential of the Sphagnum species composition to adjust to this forcing remains poorly understood. Particularly in oligotrophic fens, where the vegetation is substantially shaped by lateral hydrology (Tahvanainen, 2011; Turetsky et al., 2012), plant communities can be highly vulnerable to hydrological changes (Gunnarsson et al., 2002; Tahvanainen, 2011). Based on the validity and robustness of PMS, we believe that this model could serve as one of the first mechanistic tools to investigate the direction and rate of change in Sphagnum communities under environmental forcing. The hummock–lawn differences showed by Test 1 imply that S. magellanicum could outcompete S. fallax within a decade in a poor fen community if the water table of habitats such as lawns was lowered by 15 cm (Test 1). Although this was derived from a simplified system comprising only the two species, it highlighted the potential for rapid turnover of Sphagnum species: the hummock–lawn difference of the water table in the simulation was comparable to the expected water table drawdown in fens under a warming climate (Whittington and Price, 2006; Gong et al., 2012). The effect traits of mosses, while less studied than those of vascular plant traits, have far reaching impacts on the biogeochemistry of ecosystems such as peatlands, where mosses form the most significant plant group (Cornelissen et al., 2007). Due to the large interspecific differences of traits such as photosynthetic potential, hydraulic properties, and litter chemistry (Laiho, 2006; Straková et al., 2011; Korrensalo et al., 2017; Jassey and Signarbieux, 2019), a change in the Sphagnum community composition is likely to impact long-term peatland stability and functioning (Waddington et al., 2015). Turnover between hummock and wetter habitat species would feedback to climate, as they differ in their decomposability (Straková et al., 2012; Bengtsson et al., 2016). As hummock species produce more recalcitrant litter, the carbon bound in the system would take longer to be released back into atmosphere. In addition, the replacement of moss species that are adapted to wet conditions by hummock species is likely to result in a higher ability to maintain a carbon sink under periods of drought (Jassey, and Signarbieux, 2019).

Although efforts have been made in analytical modelling to obtain boundary conditions for equilibrium states of moss and vascular communities in peatland ecosystems (Pastor et al., 2002), the dynamic process of peatland vegetation has not been well described or included in Earth system models (ESMs). Existing ecosystem models usually consider the features of peatland moss cover as “fixed” (Sato et al., 2007; Wania et al., 2009; Euskirchen et al., 2014) or assume that they change directionally following a projected trajectory (Wu and Roulet, 2014). Chaudhary et al. (2017) presented a dynamic peatland vegetation model with a single moss PFT and four vascular PFTs so that moss productivity can vary relative to vascular plants; however, moss characteristics are fixed to a single set of values. Our modelling approach provided a way to incorporate the environmental fluctuation and the mechanisms of dynamic moss cover into peatland carbon modelling. PMS employed an individual-based approach in which each grid cell carries a unique set of trait properties; thus, shoots with favourable trait combinations in the prevailing environment are able to replace those whose trait combinations are less favourable. Moreover, the model included the spatial interactions of individuals, which can impact the sensitivity of the coexistence pattern to environmental changes (Bolker et al., 2003; Sato et al., 2007; Tatsumi et al., 2019). This mimics the stochasticity in plant responses to environmental fluctuations, which is essential for community assembly and trait filtering under environmental forcing (Clark et al., 2010). The importance of incorporating environmental fluctuations into the stochasticity of biophysiological processes is supported by tests 9 and 10. If the monthly mean climate conditions were used as input, our model failed to predict the dominance of S. magellanicum on hummocks. If the stochasticity of model parameters was omitted and only mean values were used, the model generated only single output that disregarded the randomness of environmental conditions. As these features are considered essential for “next-generation” dynamic vegetation models (DVMs; Scheiter et al., 2013), our PMS could be considered as an elemental design for future DVM development.

We see PMS as an elemental design for the future development of dynamic vegetation models for peatland ecosystems; however, there are certain uncertainties and features that should be developed further. We used a gas-exchange-based method to quantify the simultaneous changes in the capitula water potential, the water content, and the carbon uptake of Sphagnum moss capitula. It should be noted that the measurements mainly covered the changes from WC_opt towards WC_cmp (Appendix E and Fig. 3). However, the capitula water content could be higher than WC_opt at saturation (e.g. about 25–30 g g⁻¹; Schipperges and Rydin, 1998). When the RWC is high, vapour diffusion may occur mainly from the capitula surface or macropores (instead of the inside capitula). Hence, our methodology may not be suitable to reflect the water potential changes under near-saturated conditions. In our model simulations, we used the volumetric water content of the moss carpet to estimate the RWC as an approximation for wet conditions (Eq. 11). The accuracy of such an approximation for high RWC conditions remains ambiguous, and more information is still required. We assumed that tissue structure did not change during the measurement process and that the aerodynamic resistance (r_a, Eq. 3) for vapour to diffuse from the inner capitula to the headspace was constant. However, capitula drying may change leaf curvature, especially in species with slim and sparsely spread leaves (Laine et al., 2018). Such changes in the branch-leaf structure could expose more of the leaf surface to evaporation and reduce the value of r_a. Consequently, PMS could underestimate the capitula water potential towards the drying end for those species if a constant r_a is derived from the maximal evaporation rate (E_m, Eq. 5; Fig. 3c).

The water retention relationship in PMS may not sufficiently capture water potential changes at wet and dry extremes (e.g. S. magellanicum in Fig. 4c). Water retention functions developed for mineral soils (e.g. Clapp and Hornberge, 1978; van Genuchten, 1980) may not be well parameterized for peat soils and moss (nonvascular) vegetation, particularly under very dry or wet conditions. Hence, further studies are needed to improve the description of the nonlinearity of the capitula water content, as influenced by capitula morphology (e.g. capitula biomass and shoot density) and structural changes in leaves.

PMS lacks horizontal (lateral) water transport that may allow individuals from lawn species to be present on hummocks (Rydin, 1985). With additional experimental data, such as species-specific hydraulic conductivity, the current model could be improved to also quantify the horizontal water transport among neighbouring grid cells.

To conclude, PMS could successfully capture the habitat preferences of the modelled Sphagnum species. In this respect, PMS could provide fundamental support for the future development of dynamic vegetation models for peatland ecosystems. Based on our findings, capitulum water processes should be considered as a control on vegetation dynamics in future impact studies on peatlands under changing environmental conditions.

Appendix A: Calculating the community SVAT scheme (Module III)

A1 Transport of water and heat in the peat profile

Simulating the transport of water and heat in the peat profiles was based on Gong et al. (2013). Here we list the key algorithms and parameters. Ordinary differential equations governing the vertical transport of water and heat in the peat profiles were given as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{C}_{\mathrm{h}}{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathrm{h}}\left({\displaystyle \frac{\partial h}{\partial z}}+\mathrm{1}\right)\right]+S_{\mathrm{h}}\text{(A2)}& {\displaystyle }{C}_{\text{T}}{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K_{\text{T}}{\displaystyle \frac{\partial T}{\partial z}}\right)+S_{\text{T}},\end{array}$$

where t is the time step; z is the thickness of the peat layer; h is the water potential; T is the temperature; C_h and C_T are the specific capacity of water (i.e. $\partial \mathit{\theta }/\partial h$) and heat respectively; K_h and K_T are the hydraulic conductivity and thermal conductivity respectively; and S_h and S_T are the sink terms for water and energy respectively.

C_T and K_T were calculated as the volume-weighted sums from components of water, ice, and organic matter:

$$\begin{array}{}\text{(A3)}& {\displaystyle }{C}_{\text{T}}=C_{\text{water}}{\mathit{\theta }}_{\text{water}}+C_{\text{ice}}{\mathit{\theta }}_{\text{ice}}+C_{\text{org}}\left(\mathrm{1}-{\mathit{\theta }}_{\text{water}}-{\mathit{\theta }}_{\text{ice}}\right)\text{(A4)}& {\displaystyle }{K}_{\text{T}}=K_{\text{water}}{\mathit{\theta }}_{\text{water}}+K_{\text{ice}}{\mathit{\theta }}_{\text{ice}}+K_{\text{org}}\left(\mathrm{1}-{\mathit{\theta }}_{\text{water}}-{\mathit{\theta }}_{\text{ice}}\right),\end{array}$$

where C_water, C_ice, and C_org are the specific heats of water, ice, and organic matter respectively; K_water, K_ice, and K_org are the thermal conductivities of water, ice, and organic matter respectively; and θ_water and θ_ice are the volumetric contents of water and ice respectively.

For a given h, $C_{\mathrm{h}}=\partial \mathit{\theta }\left(h\right)/\partial h$ was derived from the van Genuchten water retention model (van Genuchten, 1980) as follows:

$$\begin{array}{}\text{(A5)}& \mathit{\theta }\left(h\right)={\mathit{\theta }}_{\mathrm{r}}+{\displaystyle \frac{\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}{\left[\mathrm{1}+{\left(\mathit{\alpha }\left|h^n\right|\right)}^m\right]}},\end{array}$$

where θ_s is the saturated water content; θ_r is the permanent wilting point water content; α is a scale parameter that is inversely proportional to mean pore diameter; n is a shape parameter; and $m=\mathrm{1}-\mathrm{1}/n$.

Hydraulic conductivity (K_h) in an unsaturated peat layer was calculated as a function of θ by combining the van Genuchten model with the Mualem model (Mualem, 1976):

$$\begin{array}{}\text{(A6)}& K_{\mathrm{h}}\left(\mathit{\theta }\right)=K_{\text{sat}}{S}_{\mathrm{e}}^{L_{\mathrm{e}}}\left[\mathrm{1}-{\left(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}\right)}^m\right],\end{array}$$

where K_sat is the saturated hydraulic conductivity; S_e is the saturation ratio, where $S_{\mathrm{e}}=(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}})/({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})$; and L_e is the shape parameter (L_e=0.5; Mualem, 1976).

A2 Boundary conditions and the surface energy balance

A zero-flow condition was assumed at the lower boundary of the peat column. The upper boundary condition was defined by the surface energy balance, which was driven by net radiation (R_n). The dynamics of R_n at surface x (x=0 for vascular canopy and x=1 for moss surface) was determined by the balance between incoming and outgoing radiation components:

$$\begin{array}{}\text{(A7)}& R_{{\mathrm{n}}_x}={\text{Rsn}}_{\mathrm{b},x}+{\text{Rsn}}_{\mathrm{d},x}+{\text{Rln}}_x,\end{array}$$

where Rsn_b,x and Rsn_d,x are the absorbed energy from direct and diffuse radiation respectively; and Rln_x is the absorbed net longwave radiation.

Algorithms for calculating the net radiation components were detailed in Gong et al. (2013), as modified from the methods of Chen et al. (1999). Canopy light interception was determined by the light-extinction coefficient (k_light), the leaf area index (L_c), and the solar zenith angle. The partitioning of reflected and absorbed irradiances at ground surface was regulated by the surface albedos for the shortwave (a_s) and longwave (a_l) components, and the temperature of surface x (T_x) also affects net longwave radiation:

$$\begin{array}{}\text{(A8)}& {\displaystyle }& {\displaystyle }{R}_{{\mathrm{n}}_x}={\text{Rsn}}_{\mathrm{b},x}+{\text{Rsn}}_{\mathrm{d},x}+{\text{Rln}}_x\text{(A9)}& {\displaystyle }& {\displaystyle }{\text{Rsn}}_{\mathrm{d},x}=R_{{\mathrm{s}}_{\mathrm{id},x}}\left(\mathrm{1}-a_{\mathrm{s}}\right){\text{Rln}}_x={R_{\mathrm{l}}}_{i,x}\left(\mathrm{1}-a_{\mathrm{l}}\right)-\mathit{\epsilon }\mathit{\delta }{T}_x^{\mathrm{4}},\end{array}$$

where R_sib, R_sid, and R_li are the incoming beam, diffusive, and longwave radiation respectively; ε is the emissivity ($\mathit{\epsilon }=\mathrm{1}-a_{\mathrm{l}}$); and δ is the Stefan–Boltzmann constant ($\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴).

$R_{{\mathrm{n}}_x}$ was partitioned into the latent heat flux (λE_x), the sensible heat flux (H_x), and the ground heat flux (for canopy G₁=0):

$$\begin{array}{}\text{(A10)}& {\displaystyle }{R}_{{\mathrm{n}}_x}=H_x+\mathit{\lambda }{E}_x+G_x\text{(A11)}& {\displaystyle }{G}_{\mathrm{1}}=K_{\text{T}}\left(T_x-T_{\mathrm{s}}\right)/\left(\mathrm{0.5}z\right),\end{array}$$

where T_s is the temperature of the moss carpet; and z is the thickness of the moss layer (z=5 cm).

The latent heat flux was calculated using the “interactive scheme” (Daamen and McNaughton, 2000; see also Gong et al., 2016), which is a K-theory-based, multisource model:

$$\begin{array}{}\text{(A12)}& \mathit{\lambda }{E}_x={\displaystyle \frac{\left(\mathrm{\Delta }/\mathit{\gamma }\right)A_x{r}_{\mathrm{sa},x}+\mathit{\lambda }{\text{VPD}}_{\mathrm{b}}}{r_{\mathrm{b},x}+\left(\mathrm{\Delta }/\mathit{\gamma }\right)r_{\mathrm{sa},x}}},\end{array}$$

where Δ is the slope of the saturated vapour pressure curve against air temperature; λ is the latent heat of vaporization; E is the evaporation rate; VPD_b is the vapour pressure deficit at z_b; r_b,x is the total resistance to water vapour flow, which is the sum of boundary layer resistance (r_sa,x) and surface resistance (r_ss); and A is the available energy for evapotranspiration and A_x = $R_{{\mathrm{n}}_x}$ – G_x.

The calculations of γ, λ, and VPD_b require the temperature (T_b) and vapour pressure (e_b) at the mean source height (z_b). These variables were related to the total of latent heat (∑λE_x) and sensible heat (∑H_x) from all surfaces using the Penman-type equations:

$$\begin{array}{}\text{(A13)}& {\displaystyle }\mathrm{\Sigma }\mathit{\lambda }{E}_x={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}\left(e_{\mathrm{b}}-e_{\mathrm{a}}\right)/\left(r_{\text{aero}}\mathit{\gamma }\right)\text{(A14)}& {\displaystyle }\mathrm{\Sigma }{H}_x={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}\left(T_{\mathrm{b}}-T_{\mathrm{a}}\right)/r_{\text{aero}},\end{array}$$

where ρ_aC_p is the volumetric specific heat of air; r_aero is the aerodynamic resistance between z_b and the reference height z_a and was a function of T_b accounting for the atmospheric stability (Choudhury and Monteith, 1988); and γ is the psychrometric constant ($\mathit{\gamma }={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}/\mathit{\lambda }$).

Changes in the energy balance affect the surface temperature (T_x) and vapour pressure (e_x), which further feed back to the energy availability (Eqs. A10, A12), the source-height temperature, VPD, and the resistance parameters (e.g. r_aero). The values of T_x and e_x were solved iteratively by coupling the energy balance equations (Eqs. A11–A15) with the Penman-type equations (see also Appendix B in Gong et al., 2016):

$$\begin{array}{}\text{(A15)}& {\displaystyle }\mathit{\lambda }{E}_x={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}\left(e_x-e_{\mathrm{b}}\right)/\left(r_{\mathrm{sa},x}\mathit{\gamma }\right)\text{(A16)}& {\displaystyle }{H}_x={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}\left(T_x-T_{\mathrm{b}}\right)/r_{\mathrm{sa},x},\end{array}$$

where the boundary layer resistance for ground surface (r_sa,1) and canopy (r_sa,0) were calculated following the approaches of Choudhury and Monteith (1988).

A3 Sink terms of transport functions for water and heat

The sink term S_h,i (see Eq. A11) for each soil layer i was calculated as follows:

$$\begin{array}{}\text{(A17)}& S_{\mathrm{h},i}=E_i-P_i-W_{\text{melt},i}-I_i,\end{array}$$

where E_i is the evaporative loss of water from the layer; P_i is rainfall (P_i=0 if the layer is not the topmost layer, i.e. i<1); W_melt,i is the amount of melt water added to the layer; and I_i is the net water inflow and was calibrated in Sect. 2.5.

The value of E_i was calculated as follows:

$$\begin{array}{}\text{(A18)}& E_i=f_{\text{top}}{E}_{\mathrm{0}}+f_{\text{root}}\left(i\right)E_{\mathrm{1}},\end{array}$$

where E₀ and E₁ are the evaporation rate from the ground surface and canopy (Eq. A13) respectively; f_top is the location multiplier for the topmost layer (f_top=0 in cases i>1); and f_root(i) is the fraction of fine-root biomass in layer i.

The value of W_melt,i was controlled by the freeze–thaw dynamics of the soil water and snowpack, which were related to the heat diffusion in the soil profile (Eq. A2). We set the freezing point temperature to 0 ^∘C, and the temperature of a soil layer was held constant (0 ^∘C) during freezing or melting. For the ith soil layer, the sink term (S_T) in the heat transport equation (Eq. A2) was calculated as follows:

$$\begin{array}{}\text{(A19)}& S_{\mathrm{T},i}=f_{\text{phase}}\text{max}\left(\left|T_i\right|C_{\mathrm{T},i},W_{\text{phase}}{\mathit{\lambda }}_{\text{melt}}\right),\end{array}$$

where C_T,i is the specific heat of the soil layer (Eq. A13); W_phase is the water content for freezing (W_phase=θ_w) or melting (W_phase=θ_ice); λ_melt is the latent heat of freezing; f_phase is a binarial coefficient that denotes the existence of freezing or thawing. For each time step t, we computed T_i(t) with the prior assumption that $S_{\mathrm{T},i}=$0. f_phase was then determined by whether the temperature changed across the freezing point, i.e. f_phase=1 if $T_i\left(t\right)\times T_i(t-\mathrm{1})\le \mathrm{0}$, otherwise f_phase=0.

A4 Parameterization of SVAT processes

For the calculation of the surface energy balance, we set the height and leaf area of the vascular canopy to 0.4 m and 0.1 m² m⁻² respectively, which is consistent with the scarcity of vascular canopies at the site. The aerodynamic resistance (r_aero, Eq. A14, Appendix A) for surface energy fluxes was calculated following Gong et al. (2013a). The bulk surface resistance of the community (r_ss, Eq. A13, Appendix A) was summarized from the cell-level values of r_bulk,i such that $\mathrm{1}/r_{\text{ss}}=\sum (\mathrm{1}/r_{\text{bulk},i})$. To calculate the peat hydrology and water table, peat profiles of hummock and lawn communities were set to 150 cm deep and stratified into horizontal depth layers varying from 5 cm (topmost) to 30 cm (deepest). For each peat layer, the thermal conductivity (K_T) of the fractional components, i.e. peat, water, and ice, were evaluated following Gong et al. (2013a). The bulk density of peat (ρ_bulk) was set to 0.06 g cm⁻³ below acrotelm (40 cm depth, Laine et al., 2004) and decreased linearly toward the living moss layer. The saturated hydraulic conductivity (K_sat, Eq. A6, Appendix  A) and water retention parameters (i.e. α and n, Eq. A5, Appendix A) of water retention curves were calculated as functions of ρ_bulk and the depth of the peat layer following Päivänen (1973). K_sat, α, and n for the living moss layer were adopted from the values measured by McCarter and Price (2014) for S. magellanicum carpet. The parameter values for SVAT processes are listed in Table 3.

A5 Calculation of snow dynamics

In boreal and arctic regions, the amount and timing of snow melt has a crucial impact on moisture conditions, especially in fen peatlands. Therefore, in order to have realistic spring conditions, we introduced a snowpack model, SURFEX v7.2 (Vionnet et al., 2012) into the SVAT model simulations. The snowpack model simulates snow accumulation, wind drift, compaction, and changes in metamorphism and density. These processes influenced the heat transport and freezing–melting processes (i.e. S_h and S_T, see Eqs. A1–A2 in Appendix A). In this model run, we calculate the snow dynamics on a daily basis in parallel to the SVAT simulation. Daily snowfall was converted into a snow layer and was added to ground surface. For each of the day-based snow layers (D-layers), we calculated the changes in snow density, particle morphology, and layer thicknesses. At each time step, D-layers were binned into layers with depths of 5–10 cm (S-layers) and placed on top of the peat column for SVAT model simulations. With a snow layer present, surface albedos (i.e. a_s and a_l) were modified to match those of the topmost snow layer (see Table 4 in Vionnet et al., 2007). If the total thickness of snow was less than 5 cm, all D-layers were binned into one S-layer. The thermal conductivity (K_T), specific heat (C_T), snow density, thickness, and water content of each S-layer were calculated as the mass-weighted means from the values of D-layers. Melting and refreezing tended to increase the density and K_T of a snow layer but decrease its thickness (see Eq. 18 in Vionnet et al., 2007). The fraction of melted water that exceeded the water holding capacity of a D-layer (see Eq. 19 in Vionnet et al., 2007) was removed immediately as infiltration water. If the peat layer underneath was saturated, the infiltration water was removed from the system as lateral discharge.

A6 Boundary conditions and driving variables

A zero-flow boundary was set at the bottom of the peat. At the peat surface the boundary conditions of water and energy were defined by the ground surface temperature (T₀, see Eqs. A10–A15 in Appendix A) and the net precipitation (P minus E). The profiles of layer thicknesses, ρ_bulk, and hydraulic parameters were assumed to be constant during simulation. Lateral boundary conditions were used to calculate the spreading of Sphagnum shoots among cells along the edge of the model domain so that shoots can spread across the edge of the simulation area and invade the grid cell at the boarder of the opposite side.

The model simulation was driven by the climatic variables of air temperature (T_a), precipitation (P), relative humidity (RH), wind speed (u), incoming shortwave radiation (R_s), and longwave radiation (R_l). To support the stochastic parameterization of the model and Monte Carlo simulations, the Weather Generator (Strandman et al., 1993) was used to generate randomized scenarios based on long-term weather statistics (period from 1981 to 2010) from the four closest weather stations of the Finnish Meteorological Institute. This generator had been intensively tested and applied under Finnish conditions (Kellomäki and Väisänen, 1997; Venäläinen et al., 2001; Alm et al., 2007). We also compared the simulated meteorological variables with 2 years of data measured from a Siikaneva peatland site (61^∘50 N; 24^∘10 E) that is located 10 km from our study site (Appendix C).

Appendix B: Methods and results of the empirical study on Sphagnum capitula water retention as a controlling mechanism for peatland moss community dynamics

B1 Measurement of morphological traits

To quantify morphological traits, core (size: 7 cm diameter, 50 cm² area, and a height of at least 8 cm) samples of S. fallax and S. magellanicum were collected at the end of August 2016, maintaining the natural density of the stand. Samples were stored in plastic bags in a cool room (4 ^∘C) until analysis. Eight replicates were collected for each species. For each sample, the capitulum density (D_S, shoots cm⁻²) was measured, and 10 moss shoots were randomly selected and separated into the capitula and stems (5 cm below the capitula). The capitula and stems were moistened and placed on tissue paper for 2 min to extract free-moving water, before they were weighed for the water-filled fresh weight. The samples were then dried at 60 ^∘C for at least 48 h, and the dry masses were subsequently measured. The field-water contents of the capitula (W_cf, g g⁻¹) and stems (W_sf, g g⁻¹) were then calculated as the ratio of water to dry mass for each sample. The biomass of the capitula (B_cap, g m⁻²) and living stems (B_st, g m⁻²) was calculated by multiplying the dry mass by the capitulum density (D_S). The biomass density of living stems (H_spc, g cm⁻¹ m⁻²) was calculated by dividing B_st by the length of stems.

B2 Measurement of photosynthetic traits

We measured the photosynthetic light response curves for S. fallax and S. magellanicum with fully controlled, flow-through gas-exchange fluorescence measurement systems (GFS-3000, Walz, Germany; Li-6400, LI-COR, US) under varying light levels. In 2016, measurements on field-collected samples were carried out during May and early June, which is a peak growth period for Sphagna (Korrensalo et al., 2017). Samples were collected from the field site each morning and were measured the same day at Hyytiälä field station. Samples were stored in plastic containers and moistened with peatland water to avoid changes in the plant status during the measurement. Right before the measurement, we separated Sphagnum capitula from their stems and dried them lightly using tissue paper before placing an even layer of them into a custom-made cuvette and retaining the same density as that found under natural field conditions (Korrensalo et al., 2017). The net photosynthesis rate (A, µmol g⁻¹ s⁻¹) was measured at 1500, 250, 35, and 0 µmol m⁻² s⁻¹ photosynthetic photon flux density (PPFD; Fig. 1b). The light levels were chosen based on previous investigation by Laine et al. (2011, 2015), which showed increasing A until PPFD at 1500 and no photoinhibition even at high values of 2000 µmol m⁻² s⁻¹. The samples were allowed to adjust to cuvette conditions before the first measurement and after each change in the PPFD level until the CO₂ rate had reached a steady level, otherwise the cuvette conditions were kept constant (temperature 20 ^∘C, CO₂ concentration 400 ppm, flow rate 500 µmol s⁻¹, impeller at level 5, and relative humidity of inflow air 60 %, although the relative humidity remained on average 81 % during the measurements). The time required for a full measurement cycle varied between 60 and 120 min. Each sample was weighed before and after the gas-exchange measurement and was then dried at 40 ^∘C for 48 h to determine the biomass of capitula (B_cap). For each species, five samples were measured as replicates and were made to fit a hyperbolic light-saturation curve (Larcher, 2003):

$$\begin{array}{}\text{(B1)}& A_{\mathrm{20}}=\left({\displaystyle \frac{P_{\mathrm{m}\mathrm{20}}\times \text{PPFD}}{{\mathit{\alpha }}_{\text{PPFD}}+\text{PPFD}-R_{{\mathrm{s}}_{\mathrm{20}}}}}\right)\times B_{\text{cap}},\end{array}$$

where subscript 20 denotes the variable value measured at 20 ^∘C; R_s is the mass-based dark respiration rate (µmol g⁻¹ s⁻¹); P_m is the mass-based rate of maximal gross photosynthesis (µmol g⁻¹ s⁻¹); and α_PPFD is the half-saturation point (µmol m⁻² s⁻¹), i.e. the PPFD level where half of P_m is reached. The measured morphological and photosynthetic traits are listed in Table 2.

Figure B1 Measured light response curves for S. magellanicum and S. fallax.

B3 Drying experiment

To link the water retention and photosynthesis of Sphagnum capitula, we performed a drying experiment using a GFS-3000 system to measure covariations of capitulum water potential (h, cm water), water content (W_cap, g g⁻¹), and A (µmol g⁻¹ s⁻¹). For both species, four mesocosms were collected in August 2018 and transported to the laboratory at UEF Joensuu, Finland. Capitula were harvested and wetted by water from the mesocosms. The capitula were then placed gently onto a piece of tissue paper for 2 min before being placed into the same cuvette as that used in the previous photosynthesis measurement. The cuvette was then placed into the GFS-3000 and measured under constant PPFD (1500 µmol m⁻² s⁻¹), temperature (293.2 K), inflow air (700 µmol s⁻¹), CO₂ concentration (400 ppm), and relative humidity (40 %) conditions. Measurement was stopped when A dropped below 10 % of its maximum. Each measurement lasted between 120 and 180 min. Each sample was weighed before and after the gas-exchange measurement and then dried at 40^∘C for 48 h to determine the biomass of capitula (B_cap).

The GFS-3000 records the vapour pressure (e_a, kPa) and the evaporation rate (E, g s⁻¹) simultaneously with A every second (Walz, Germany, 2012). The changes in W_cap with time (t) were calculated as follows:

$$\begin{array}{}\text{(B2)}& \text{RWC}\left(t\right)=\left(W_{\text{pre}}-B_{\text{c}}-{\sum }_{t=\mathrm{0}}^{t}E\left(t\right)\right)/B_{\text{c}}\end{array}$$

We assumed that the vapour pressure at the surface of water-filled cells equaled the saturation vapour pressure (e_s), and the vapour pressure in the headspace of the cuvette equaled that in the outflow (e_a). Thus, the vapour pressure in capitula pores (e_i) can be calculated based on the following gradient-transport function (Fig. B2a):

$$\begin{array}{}\text{(B3)}& \mathit{\lambda }E\left(t\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}}{\mathit{\gamma }}}{\displaystyle \frac{\left(e_i\left(t\right)-e_{\mathrm{a}}\left(t\right)\right)}{r_{\mathrm{a}}\left(t\right)}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}}{\mathit{\gamma }}}{\displaystyle \frac{\left(e_{\mathrm{s}}-e_i\left(t\right)\right)}{r_{\mathrm{s}}\left(t\right)}},\end{array}$$

where λ is the latent heat of vaporization; γ is the slope of the saturation vapour pressure–temperature relationship; r_a is the aerodynamic resistance (m s⁻¹) for vapour transport from the inter-leaf volume to the headspace; and r_s is the surface resistance of vapour transport from the wet leaf surface to the inter-leaf volume. Thus, the bulk resistance for evaporation (r_bulk) was calculated as r_a+r_s.

We assumed that the structures of the tissues and pores did not change during the drying process and that r_a was constant during each measurement. A tended to increase with time t until it peaked (A_m) and then subsequently decreased (Fig. 2b). The point A=A_m implied the water content where further evaporative loss would start to drain the cytoplasmic water, leading to a decrease in A. The response of A to W_cap was fitted as a second-order polynomial function (Robroek et al., 2009) using data from $t_{A_{\mathrm{m}}}$ to t_n:

$$\begin{array}{}\text{(B4)}& f_A\left(W_{\text{cap}}\right)=a_{\text{W0}}+a_{\text{W1}}\times W_{\text{cap}}+a_{\text{W2}}\times W_{\text{cap}}^{\mathrm{2}},\end{array}$$

where a_W0, a_W1, and a_W2 are parameters, and $f_A\left(W_{\text{cap}}\right)=A/A_{\mathrm{m}}$. For each replicate, the optimal water content for photosynthesis (W_opt) was derived from the peak of the fitted curve (Eq. 4). The capitulum water content at the compensation point W_cmp, where the rates of gross photosynthesis and respiration are equal, can be calculated from the point A=0.

Figure B2 Conceptual schemes of (a) the cuvette setting and resistances, (b) the covariations of net photosynthesis and W_cap, and (c) the covariations of evaporation and vapour pressure in the headspace during a measurement. The symbols utilized in the figure are as follows: e_a is the vapour pressure in the headspace of the cuvette (kPa); e_i is the vapour pressure in the branch-leaf structure of the capitula; e_s is the vapour pressure at the surface of wet tissues; r_a is the aerodynamic resistance of vapour diffusion from the inner capitula to the headspace; r_s is the surface resistance of vapour diffusion from the wet tissue surface to the inner capitula space; A is the net photosynthesis rate (µmol m⁻² s⁻¹); A_m is the maximal net photosynthesis rate (µmol m⁻² s⁻¹); W_cap is the water content of the capitula (g g⁻¹); W_opt is W_cap at A=A_m; W_cmp is W_cap at A=0; E is the evaporation rate (mm s⁻¹).

Similarly, the evaporation rate (E) increased from the start of the measurement until the maximum evaporation E_m and then decreased (Fig. B2c). The point E=E_m implied the time when the wet capitulum tissues underwent maximum exposure to the air flow. Therefore, r_a was estimated as the minimum of bulk resistance using Eq. (B5) and assuming e_i(t)≈e_s when E(t)=E_m:

$$\begin{array}{}\text{(B5)}& r_{\mathrm{a}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}}{\mathit{\gamma }}}{\displaystyle \frac{\left(e_{\mathrm{s}}-e_{\mathrm{a}}\left(t\right)\right)}{\mathit{\lambda }{E}_{\mathrm{m}}}}\end{array}$$

Based on the calculated e_i(t), we were able to derive the capitulum water potential (h) following the equilibrium vapour pressure method (e.g. Price et al, 2008; Goetz and Price, 2015):

$$\begin{array}{}\text{(B6)}& h={\displaystyle \frac{RT}{M_g}}\ln \left({\displaystyle \frac{e_i}{e_{\mathrm{s}}}}\right)+h_{\mathrm{0}},\end{array}$$

where R is the universal gas constant (8.314 J mol⁻¹ K⁻¹); M the molar mass of water (0.018 kg mol⁻¹); g is the gravitational acceleration (9.8 N kg⁻¹); e_i∕e_s is the relative humidity; and h₀ is the water potential due to the emptying of free-moving water before measurement (set to 10 kPa according to Hayward and Clymo, 1982).

Table B1 Species-specific morphological, photosynthetic, and water retention traits from S. magellanicum and S. fallax. Trait values (mean ± standard deviation) and ANOVA statistics, F and p values, are given to compare the mean values of traits for the two species.

^a Soil core measurement, sample n=5. ^b Cuvette gas-exchange measurement, sample n=4. ^∗ The difference between the means is significant (P<0.05). ${}^{\ast \ast }$ The difference between the means is very significant (P<0.01).

Table B2 Parameter estimates of the linear model for the log₁₀-transformed capitulum water potential (h) for S. fallax and S. magellanicum. The estimated value, the standard error (SE), and the test statistics (p values) are given for the predictors of the models. The predictors are the capitulum biomass (B_cap), the capitulum density (D_S), the capitulum water content (W_cap), the interaction of the capitulum biomass and water potential (B_cap×W_cap), the interaction of the capitulum biomass and capitulum density (D_S×W_cap), the interaction of the capitulum density and water potential (D_S×W_cap), and the interaction of the capitulum biomass, capitulum density, and water potential ($B_{\text{cap}}\times D_{\text{S}}\times W_{\text{cap}}$). All coefficient values are significantly different from zero (p<0.001).

Table B3 Parameter estimates of the linear model for the log₁₀-transformed capitulum evaporative resistance (r_bulk) for S. fallax and S. magellanicum. The estimated value, the standard error (SE), and the test statistics (p values) are given for the predictors of the models. The predictors are the capitulum biomass (B_cap), the capitulum density (D_S), the water potential (h), the interaction of capitulum biomass and water potential (B_cap×h), the interaction of capitulum biomass and capitulum density (D_S×h), the interaction of capitulum density and water potential (D_S×h), and the interaction of capitulum biomass, capitulum density, and water potential ($B_{\text{cap}}\times D_{\text{S}}\times h$). All coefficient values are significantly different from zero (p<0.001).

B4 Statistical analysis

The light response curve (Eq. B1) and the response function of A∕A_m to W_cap changes (Eq. B4) were fitted using the nlme package in R (version 3.1). The obtained values of the shape parameters a_W0, a_W1, and a_W2 (Eq. 4) were then used to calculate W_opt ($W_{\text{opt}}=-\mathrm{0.5}{a}_{\text{W1}}/a_{\text{W2}}$) and W_cmp $(W_{\text{cmp}}=\mathrm{0.5}[-a_{\text{W1}}-(a_{\text{W1}}^{\mathrm{2}}-\mathrm{4}{a}_{\text{W0}}{a}_{\text{W2}}{)}^{\mathrm{0.5}}]/a_{\text{W2}})$. We subsequently applied an ANOVA to compare S. magellanicum with S. fallax with respect to the traits obtained from the field sampling (i.e. structural properties such as B_cap, D_S, H_spc, W_cf, and W_sf) and from the gas-exchange measurements (i.e. P_m20, R_s20, W_opt, W_cmp, and r_bulk), using R (version 3.1).

The measured values of the capitulum water potential (h) were log₁₀ transformed and related to the variations in W_cap, B_cap, and D_S using a linear model. Similarly, a linear model was established to quantify the response of bulk resistance for evaporation (r_bulk) (log₁₀ transformed) to the variations in h, B_cap, and D_S. The linear regressions were based on statsmodels (version 0.9.0) in Python (version 2.7), as supported by the NumPy (version 1.12.0) and pandas (version 0.23.4) packages.

Figure B3 Responses of (a) the capitulum water potential, (b) the bulk resistance of evaporation, and (c) the net photosynthesis to changes in the capitulum water content (W_cap) of two Sphagnum species typical to hummocks (S. magellanicum, black) and lawns (S. fallax, red). As the measured results are based on the drying experiment starting with fully wetted capitula characteristic for both species, the x axis is shown from high to low W_cap. The values predicted in panels (b) and (c) are based on linear models with the parameter values listed in Tables B2 and B3 and predictor values from the drying experiment.

B5 Results of the empirical measurements

The two Sphagnum species differed with respect their structural properties (Table B1). The lawn species S. fallax had a looser structure than the hummock species S. magellanicum, as seen by the lower capitulum density (D_S) and specific height (H_spc) in S. fallax than in S. magellanicum (P<0.05, Table. B1). Moreover, under the conditions prevailing at the study site, S. fallax mosses were dryer than S. magellanicum; the field-water contents of the S. fallax capitulum (W_cf) and stem (W_sf) were 40 % and 46 % lower than S. magellanicum (P<0.01, Table B1) respectively. The different density of the capitulum of the two species, which also differed in their capitulum size, led to similar capitulum biomass (B_cap) (P=0.682) between S. fallax with small capitulum and S. magellanicum with large capitulum. Unlike the structural properties, maximal CO₂ exchange rates (P_m20 and R_s20) did not differ between the two species (Table B1).

The drying experiment demonstrated how the capitulum water content regulated capitulum processes in both of the Sphagnum species studied (Fig. B3). Decreasing capitulum water content (W_cap) led to a decrease in the water potential (h), and the responses of h to W_cap varied among replicates (Fig. 3a). The values of W_cap for S. fallax were generally lower than those for S. magellanicum under the same water potentials. The fitted linear models explained over 95 % of the variation in the measured h for both species (Table B2), although the fitted responses of h to W_cap were slightly smoother than the measured responses, particularly for S. magellanicum (Fig. B3a). The responses of h to W_cap were significantly affected by the capitulum density (D_S), the capitulum biomass (B_cap), and their interactions with W_cap (Table B2).

Decreasing capitulum water content (W_cap) and water potential (h) were associated with increasing bulk resistance for evaporation (r_bulk, Fig. B3b), although the sensitivity of r_bulk to h changes varied by replicate. The values of r_bulk from S. fallax were largely lower than those from S. magellanicum when the capitulum water contents of the two species were similar. The fitted linear models explained the observed variations in the measured r_bulk well for both species (Fig. 2b, Table B3). The variation in the response of r_bulkto h was significantly affected by the capitulum density (D_S), the capitulum biomass (B_cap), and their interactions with h (Table B3).

Decreasing capitulum water content (W_cap) slowed the net photosynthesis rate down (Fig. B2c), as represented by the decreasing ratio of A∕A_m. S. fallax required a lower capitulum water content (W_cap) than S. magellanicum to reach photosynthetic maximum and the photosynthetic compensation point. However, the ranges of the capitulum water content from the photosynthetic maximum (W_opt) or field capacity (W_fc) to that at the compensation point (W_cmp) were smaller for S. fallax than for S. magellanicum. Hence, S. fallax had narrower transition zone for photosynthesis to respond to drying compared with S. magellanicum.

Appendix C: Comparisons of meteorological variables simulated by the Weather Generator and those measured from a Siikaneva peatland site (Integrated Carbon Observation System site located 10 km from the Lakkasuo study site)

Figure C1 Comparisons of meteorological variables simulated by the Weather Generator and those measured from a Siikaneva peatland site. The variables include (a) cumulative precipitation (mm), (b) incoming shortwave radiation (W m⁻²), (c) air temperature (^∘C), and (d) relative humidity (%). These variables were measured and simulated at a 30 min timescale. The measurements were carried out from 2012 to 2013. Details about the site and measurements are given in Alekseychik et al. (2018). The measured seasonal dynamics of the meteorological variables were generally in line with the 95 % confidence intervals (CI) of the simulated values, which were calculated based on Monte Carlo simulations (n=5).

Appendix D: Comparisons of the seasonal water table measured from the study site and the values simulated based on calibrated net inflow

Figure D1 Comparison of the seasonal water table (WT) measured at the Lakkasuo study site and the values simulated by the calibrated PMS. WT values were sampled weekly from the lawn habitats both in the field and in the model output. The weekly mean WT was measured during 2001, 2002, 2004, and 2016. The modelled means and standard deviations (SD) of the WT were based on 20 Monte Carlo simulations. The simulated seasonality of the mean WT generally followed the measured trends. The calibration reduced the sum of squared error (SSE, Eq. 12) from 199.5 ($a_{\text{N}}=b_{\text{N}}=\mathrm{0}$) to 117.3. The calibrated values for a_N and b_N were $-\mathrm{5.3575}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{4.7599}\times {\mathrm{10}}^{-\mathrm{5}}$ respectively (Eq. A18).

Appendix E: List of symbols and abbreviations.

Symbol	Description	Unit
A	Net photosynthesis rate	(µmol m−2 s−1)
Am	Maximal net photosynthesis rate	(µmol m−2 s−1)
αimm	Temperature constant for NSC immobilization	(–)
αPPFD	Half-saturation point of PPFD for photosynthesis	(µmol m−2 s−1)
Bcap	Capitulum biomass	(g m−2)
CT	Specific heat	(J K−1 kg−1)
DS	Capitulum density	(shoots cm−2)
dH	Annual vertical growth of Sphagnum mosses	(cm)
dWT	Hummock–lawn differences in the water table	(cm)
E	Rate of evaporation	(cm time step−1)
fW	Water content multiplier on the photosynthesis rate	(–)
fT	Temperature multiplier on the photosynthesis rate	(–)
h	Water potential	(cm)
Hc	Shoot height of Sphagnum mosses	(cm)
Hcap	Height of capitula	(cm)
Hspc	Biomass density of living Sphagnum stems	(g m−2 cm−1)
I	Rate of net inflow water	(cm)
kimm	Specific immobilization rate	(g g−1)
JDthaw	Julian day when thawing completed	(–)
Kh	Hydraulic conductivity of the peat layer	(cm s−1)
Km	Hydraulic conductivity of the moss layer	(cm s−1)
Ksat	Saturated hydraulic conductivity	(cm s−1)
KT	Thermal conductivity	(W m−1 K−1)
lc	Width of a grid cell in the simulation	(cm)
MB	Immobilized NSC to biomass production	(g)
NSCmax	Maximal NSC concentration in Sphagnum biomass	(g g−1)
P	Precipitation	(cm)
Pm	Mass-based rate of maximal gross photosynthesis	(µmol g−1 s−1)
PPFD	Photosynthetic photon flux density	(µmol m−2 s−1)
ρbulk	Bulk density of peat	(g cm−3)
raero	Aerodynamic resistance	(s m−1)
rbulk	Cell-level bulk surface resistance	(s m−1)
rss	Bulk surface resistance of community	(s m−1)
RH	Relative humidity	(%)
Rs	Mass-based respiration rate	(µmol g−1 s−1)
Rs	Incoming shortwave radiation	(W m−2)
Rl	Incoming longwave radiation	(W m−2)
Sc	Area of a cell in the model simulation	(m2)
simm	Multiplier for the temperature threshold	(–)
Svi	Evaporative area of a cell i	(cm2)
T	Capitulum temperature	(∘C)
Ta	Air temperature	(∘C)
Topt	Reference temperature of respiration (20 ∘C)	(∘C)
u	Wind speed	(m s−1)
Wcap	Capitulum water content	(g g−1)
Wcmp	Capitulum water content at the compensation point	(g g−1)
Wmax	Maximum water content of capitula	(g g−1)
Wopt	Optimal capitulum water content for photosynthesis	(g g−1)
Wcf	Field-water contents of the Sphagnum capitulum	(g g−1)
Wsf	Field-water contents of the Sphagnum stem	(g g−1)
WTm	Measured multiyear mean of the weekly water table	(cm)
WTs	Simulated multiyear mean of the weekly water table	(cm)
zm	Thickness of the living moss layer	(cm)
θm	Volumetric water content of the moss layer	(–)
θr	Permanent wilting point water content	(–)
θs	Saturated water content	(–)

Abbreviations
Γ	Learning rate of gradient decedent algorithms
D-layer	Daily based snow layer
ICOS	Integrated Carbon Observation System
JD	Julian day
NSC	Nonstructural carbon
PMS	Peatland Moss Simulator
RWC	Capitulum water content
SD	Standard deviation
SE	Standard error
SSE	Sum of squared error
SVAT	Soil–vegetation–atmosphere transport
WT	Water table