Biotic and Abiotic Drivers of Peatland Growth and Microtopography: A Model Demonstration

Peatlands are important carbon reserves in terrestrial ecosystems. The microtopography of a peatland area has a strong influence on its carbon balance, determining carbon fluxes at a range of spatial scales. These patterned surfaces are very sensitive to changing climatic conditions. There are open research questions concerning the stability, behaviour and transformation of these microstructures, and the implications of these changes for the long-term accumulation of organic matter in peatlands. A simple two-dimensional peat microtopographical model was developed, which accounts for the effects of microtopographical variations and a dynamic water table on competitive interactions between peat-forming plants. In a case study of a subarctic mire in northern Sweden, we examined the consequences of such interactions on peat accumulation patterns and the transformation of microtopographical structure. The simulations demonstrate plausible interactions between peatland growth, water table position and microtopography, consistent with many observational studies, including an observed peat age profile from the study area. Our model also suggests that peatlands could exhibit alternative compositional and structural dynamics depending on the initial topographical and climatic conditions, and plant characteristics. Our model approach represents a step towards improved representation of peatland vegetation dynamics and net carbon balance in Earth system models, allowing their potentially important implications for regional and global carbon balances and biogeochemical and biophysical feedbacks to the atmosphere to be explored and quantified.

variations and a dynamic water table on competitive interactions between peat-forming plants. In a case study of a subarctic mire in northern Sweden, we examined the consequences of such interactions on peat accumulation patterns and the transformation of microtopographical structure. The simulations demonstrate plausible interactions between peatland growth, water table position and microtopography, consistent with many observational studies, including an observed peat age profile from the study area. Our model also suggests that peatlands could exhibit alternative compositional and structural dynamics depending on the initial topographical and climatic conditions, and plant characteristics. Our model approach represents a step towards improved representation of peatland vegetation dynamics and net carbon balance in Earth system models, allowing their potentially important implications for regional and global carbon balances and biogeochemical and biophysical feedbacks to the atmosphere to be explored and quantified.

INTRODUCTION
Northern peatlands are one of the biggest carbon reserves among terrestrial ecosystems (Yu and others 2010). They have low net primary productivity (NPP) relative to other ecosystems globally, but due to an even lower decomposition rate they have sequestered considerable amounts of carbon over the course of centuries (Clymo 1991;Thormann and Bayley 1997;Frolking and others 2001). Northern peatlands have accumulated carbon at a rate of 15-30 g C m -2 y -1 in the last 5000-10,000 years (Yu and others 2009;Loisel and others 2014). Overall, the imbalance between the rate of NPP and microbial decomposition rate has created a repository of approximately 200-550 PgC covering about 3.5 million km 2 of peatland areas (Gorham 1991;Turunen and others 2002;Yu 2012).
Global climate models have projected a warming of more than 5°C by the end of the century in northern peatland areas (Christensen and others 2007;IPCC 2013). The magnitude of this warming will be exacerbated by the associated arctic climatecarbon feedbacks (Bridgham and others 1995;Davidson and Janssens 2006;Zhang and others 2014). Constraints on biological activity imposed by low temperature would be reduced, accelerating plant productivity as well as decomposition rates (Klein and others 2013). Potentially, the resulting shift in the balance between production and decomposition might be sufficient to influence the existing sink capacity of these peatlands (Wieder 2001;Ise and others 2008;Fan and others 2013). Studies have also highlighted the implications of alterations in precipitation patterns and rapid rates of permafrost degradation on the peatland carbon cycle and climate (Christensen and others 2004;Å kerman and Johansson 2008;Bragazza and others 2013). These altered patterns have the potential to modify local and regional hydrology and moisture balances, which could in turn affect the vegetation composition and carbon balance of many northern peatlands. For instance, abrupt changes in environmental conditions may create unfavourable settings for existing plant species and provide opportunities for new (non-resident) species to flourish in those areas (Malmer and others 2005;Johansson and others 2006;Zhang and others 2013).
Hummock and hollow microtopography is a distinctive feature of many northern peatland ecosystems, with microtopographical position playing a critical role in carbon balance (Gorham 1991;Weltzin and others 2001;Pouliot and others 2011). Local climate conditions, surface hydrology and vegetation cover are the main factors controlling the dynamics of patterned surfaces. These microformations develop distinctive conditions with plant species, hydrology, nutrient status, plant productivity and decomposition rates varying systematically at a range of spatial scales (Bridgham and others 1995;Waddington and Roulet 1996;Nungesser 2003;Rietkerk and others 2004;Eppinga and others 2008). Changes in regional climatic conditions could have a profound impact on these microformations, modifying the peatland carbon balance from microscale (1-10 m) to macroscale (> 10 4 m). For example, in a well-studied peatland in subarctic northern Sweden, Stordalen mire, permafrost underlying elevated areas is being degraded as a result of recent climate warming, with an increase in wet depressions modifying the vegetation composition and overall carbon sink capacity of the mire (Christensen and others 2004;Malmer and others 2005;Johansson and others 2006).
The microtopography of peatlands comprises hummocks, hollows and intermediate areas. Hollows are dominated by tall productive graminoids, while dwarf shrubs are favoured in hummock areas where the water table position (WTP) remains relatively low. Between these two extremes, different species of moss can thrive in distinctive hydrological gradients-dry, wet and intermediate (where the WTP remains close to the surface) (Seppa 2002;Malmer and others 2005;Johansson and others 2006;Trudeau and others 2014). These plant species accumulate carbon at different rates according to their morphological and structural make-up (Aerts and others 1999). Early studies such as those of VonPost and Sernander (1910) and Osvald (1923) postulated that deposition of organic plant material at dissimilar rates leads to a cyclic transformation of hummocks into hollows and vice versa. Later studies have challenged this ''cyclical succession'' hypothesis, emphasizing instead the importance of variable decomposition rates on the microtopography of peatlands (Tolonen 1971;Barber 1981;Seppa 2002). Litter derived from plant species growing in hollow areas decomposes more quickly than that derived from the vegetation of hummocks and intermediate areas, with the result that peatland microtopography can remain stable with little alterations for long periods of time, provided that climate does not change over that period. By contrast, studies in some peatland areas have found that spatial heterogeneity or surface roughness decreases over time due to differential rates of peat accumulation in hummocks and hollows, leading to a flatter surface (Clymo and Hayward 1982;Hayward and Clymo 1983). These varied and contradictory observations on growth patterns and surface microstructure dynamics complicate the construction of suitable models of the responses of peatland microtopography, vegetation cover and carbon dynamics at the regional scale, for instance in climate change impact (McGuire and others 2012) and Earth system modelling (Zhang and others 2013) studies.
In this study, we propose and demonstrate a novel two-dimensional (2-D) model representation of peat microsite dynamics and will address the issues of whether the microtopography is static or successional and whether the total surface roughness changes over time. Earlier, similar models have generally focused on peat accumulation and decomposition patterns, overlooking the tight coupling between vegetation dynamics and microtopographical transformation and their role in peatland carbon balance (Frolking and others 2001;Bauer 2004;Wania et al. 2009aWania et al. , 2009bFrolking and others 2010). The proposed model fills this gap by testing the adequacy of simple but empirically defensible assumptions about the dynamic responses of the major peatland plant functional groups to microtopography and water table variations, taking into account many of the major feedbacks of plant compositional shift on hydrology, decomposition and peat accumulation. The model is designed to allow it to be embedded within the biogeochemistry component of a regional Earth system (climate) model, accounting for the effects of vegetation dynamics on peatland carbon balance and vegetation cover at much larger (regional) scales (Smith and others 2011). With this goal in mind, we choose to adopt a parsimonious approach by restricting the number of inputs required by the model (common climate variables and primary productivity estimates) and adopting process representations that do not rely on the availability of detailed local-scale measurements (for example, soil temperature calculations). We employ the model to address the relative merits of alternative theories as to the origin and stability of peatland microtopography, in turn shedding light on the potential influence of the explored mechanisms on peatlands and their role in the global carbon balance under a changing climate. We also perform sensitivity analyses to determine the impact of model parameters and initialization on simulated microtopography, vegetation and longterm peatland carbon balance.

Model Overview
The structure of the model is presented schematically in Figure 1. Relative abundances of shrubs, graminoids and mosses in adjacent patches in a shared landscape are represented, changing over time in response to small-scale hydrology, productivity and decomposition rates that govern the rate of peat accumulation (or depletion) in each patch. These processes are assumed to determine the growth pattern, behaviour and transformation of microtopographical structures. The key formulations are adopted from previously developed peat growth models (see below). The novel feature is the inclusion of dynamic vegetation composition and the resulting impacts on peat accumulation or loss and small-scale hydrology. The input variables are annual NPP at the mean landscape level (distributed based on plant productivity), daily temperature and precipitation. The foci of the study are the peat deposition and transformations of surface microstructures during the last 1000 years. Due to this brief period of peat accumulation history considered (Stordalen peat inception has been estimated to be 4700 cal. BP-Kokfelt and others (2010)), little variation in peat bulk density with respect to depth is expected, so a constant dry bulk density is assumed in this study (Table 1).

Peat Accumulation and Decomposition
Annual peat accumulation and decomposition follow Clymo (1984), others (2001), andBauer (2004). The acrotelm and catotelm are two functionally distinct layers typical of most peatland soils. The acrotelm is made up of the comparatively aerated upper layers in the peat column, which plays the main role in determining vegetation composition. Vegetation growth results in deposits of organic material (litter) being transferred into the acrotelm. The water table fluctuates in the acrotelm, depending on microtopography, rainfall, snowmelt, evapotranspiration and run-off, with some areas remaining drier while other parts remaining completely saturated. Due to this uneven wetness, litter decomposes aerobically as well as anaerobically in the acrotelm (Clymo 1991;Frolking and others 2002). The catotelm exists below the permanent annual WTP and remains waterlogged throughout the year, creating anoxic conditions which in turn attenuate the decomposition rate and promote peat accumulation.
This model implicitly divides the total peat column into two parts-acrotelm and catotelm-demarcated by WTP. (WTP takes negative (positive) values when the WTP is below (above) the peat surface.) Annually, plant biomass is deposited as a new layer over previously accumulated peat layers in each patch (see below) based on the fractional projective cover (FPC) of the vegetation across the assumed landscape. Decomposition transforms the litter biomass in each layer into peat as time passes (Clymo 1991). The rate of change in peat mass is the total peat production minus total peat loss due to decomposition, modelled as: where M is the total peat mass (kg C m -2 ), A is the total peat input (kg C m -2 yr -1 ), and K (y -1 ) is the decomposition rate (see equation 3). Peat height may be derived from M using the assumed constant value for bulk density shown in Table 1. Peat decomposition is simulated on a daily time step based on decomposability which varies depending on the source plant type (Table 1). Litter generated from mosses (all types) has relatively low decomposition rates, shrub litter decomposes at a relatively higher rate, and graminoid litter has the highest decomposition rate; initial decomposition rates (k o -see equation 2) of 0.1, 0.05 and 0.075 y -1 for wet, intermediate and dry mosses species, respectively, and 0.15 and 0.25 y -1 for shrubs and graminoids, respectively, were assumed (Aerts and others 1999;Frolking and others 2001). The turnover rates for different source plant types decline and converge over time, with the rate of dampening being calculated using a simplified reduction equation proposed by Clymo and others (1998): where k o is the initial decomposition rate, the decay parameter n is the decay coefficient, m o is the initial mass, and m t is the mass remaining at some point in time (t). Peat water content (h) and soil temperature (T s ) have multiplicative effects on the daily decomposition rate (K i ) in each layer following Ise and others (2008) and Lloyd and Taylor (1994): where T m k i is the decomposition rate of layer i given in equation 2, and T m and W m are the temperature and moisture multipliers, respectively. We assumed that at field capacity (h opt ), peat decomposes very quickly, but when the area is waterlogged its decomposition rate decreases. We allowed the peat to decompose in very dry conditions when the WTP drops below -40 cm and water content goes below 0.01 in the peat layers.
; h ! 0:01 and h h opt b;h<0:01 and WTP < À 40 where h is the volumetric peat water content; h opt is the field capacity (0.6) and optimum volumetric water content when W m becomes 1.0; a is a parameter that affects the shape of the dependency of decay, set to 5; and b (0.6) is a minimum decomposition rate during very dry conditions when WTP goes below -40 cm. The decomposition rate is exponentially affected by soil temperature: where T i (°C) is the peat temperature in peat layer (i) and Z is a parameter affecting the slope of the exponential function (Table 1)

Hydrology
The hydrology module simulates the daily WTP which in turn determines the local vegetation cover and affects decomposition rates. A traditional water bucket scheme is adopted: where W is the total water input, P is the precipitation, ET is the evapotranspiration rate, ROF is the surface run-off, and LF is the lateral flow within the landscape depending upon the relative position of the patch. Using this scheme, peat water content (h) is updated daily in each layer of each patch (see below). Precipitation is the major source of water in ombrotrophic bogs and provides the water input to the system. Precipitation comes in the form of rain or snow depending upon the daily surface air temperature (T). When the temperature falls below the freezing point (0°C assumed), water is stored in the snow pack above the peat layers. Snow melts when the temperature rises above the freezing point, and melt rate is also influenced by daily precipitation (Choudhury and others 1998): where M is the daily total snow melt (mm), P is the daily precipitation (mm), T is the daily surface temperature (°C), T s = 0.0 is the maximum temperature (°C) at which precipitation remains snow, and SP is the current snowpack (mm). Evapotranspiration and run-off remove water from the peat column. In Stordalen, summer evaporation can reach 2-3 mm day -1 while remaining around 0.5-1 mm day -1 in winter (Rosswall and others 1975). Evapotranspiration is an empirical decreasing function of WTP (Frolking and others 2010) and is calculated using: where E o is the maximum evapotranspiration rate (Table 1) and WTP is expressed in mm. Run-off is computed following a relationship modified from Wania and others (2009a): maximum run-off can reach up to 2-3 mm day -1 when the WTP reaches + 15 cm above the surface, and there is negligible water loss from run-off below a WTP of -500 mm.
Patches can hold water up to + 15 cm above the surface giving suitable conditions for graminoids to flourish (see below). Warmer and drier conditions deepen the WTP, which in turn increases acrotelm depth. This prolongs peat residence time in the acrotelm, reducing the amount of carbon that passes down to the catotelm. Peat layers above the WTP in each patch are assumed to remain unsaturated.

Patch Structure and Lateral Water Flow
The model is run for a landscape consisting of 50 random patches of uneven height representing small-scale spatial heterogeneity typically present in peatlands. Each patch has its unique vegetation cover, hydrology, productivity and decomposition rates. Soil water across the landscape (that is, all patches) is aggregated and redistributed among patches each day through lateral flow effecting in a simple way horizontal water flow between patches. The elevated areas lose water based on the mean landscape-level WTP, whereas depressed areas receive water from the adjacent patches. In normal conditions, elevated areas have deeper WTP, whereas hollows remain waterlogged which in turn determines the vegetation, productivity and decomposition dynamics in each patch. We calculate the landscape WTP and add and remove the amount of water from each patch required to match the landscape WTP: where MWTP is the mean WTP across all the patches, PWTP i is the water table position in individual patches (i), and n is the total number of patches. The water to be added to or removed from each patch with respect to mean water table position (MWTP) in each patch, that is, lateral flow (LF), is given by: where DWTP i is the difference in the patch (i) and MWTP and LF i is the total water to be added or removed with respect to MWTP in each patch (i). If the WTP is below the surface then the total water transfer is calculated by the difference in WTP (water heights) multiplied by average porosity (U a ), whereas when the WTP is above the surface then U a is not included in the calculation. The daily water balance is implemented before the exchange of water between patches.

Vegetation Dynamics
Vegetation dynamics are affected both by autogenic (WTP and carbon mineralization) and by allogenic (temperature and precipitation) factors influencing plant litter quality and quantity, microtopography and the total amount of peat. In this study, a fractional area approach is used to determine and update annually the percentage of patch area each plant type occupies based on annual WTP. Five plant types-shrubs (Sh), graminoids (Gr), hummocks, hollows and lawn mosses (Mhu, Mho and Mim, respectively)-are considered in this study. . We assume a sigmoid relationship between annual WTP and plant productivity reduction (Yin and others 2003), expressed here through the use of a multiplicative FPC reduction term, applied annually to the FPC values of the plant types ( Figure 2). Annual WTP influences productivity differentially for each of the plant types. Mho experiences optimal productivity (no FPC reduction) when annual WTP is between + 10 cm above and -5 cm below the surface, declining at higher and lower WTPs, whereas Mim and Mhu optimal productivity is between -5 to -20 cm and -20 to -35 cm, respectively. They approach extinction (FPC reduction = 0.99) when annual WTP exceeds 10 cm above or below the optimal condition. Graminoids are characteristic of saturated conditions, experiencing optimal productivity at annual WTP > + 10 cm and approaching extinction at a WTP of -5 cm or more. At -35 cm WTP shrubs begin to occur, coming to dominate the entire patch when WTP remains below -50 cm ( Figure 2): where FPC x is the FPC for plant type x: graminoid, shrub or mosses, AWTP is the annual mean WTP, and a, b and c are parameters defined in Table 2 for  each plant type. There is, however, a lag in the system, whereby WTP changes cause marginal shifts in the relative FPCs of the plant types in each patch, each year, through multiplication of the FPC reduction term in Figure 2. The FPC of any plant type is constrained never to drop below 10 -5 allowing for recovery under suitable conditions. Emulating lags due to population dynamics and community interactions in real peatland ecosystems, the model does not determine a new equilibrium vegetation composition instantaneously. This mechanism allows vegetation to withstand short-term climatic fluctuations and also to maintain a smooth transition between the composition of plant species, expected as an outcome of the lagged effects of competition on the relative growth and mortality of the five plant types. Biomass for each plant type is initialized randomly in the first year of a simulation. FPC summed across the five plant types is constrained to 1, representing full vegetation coverage with no bare ground. Plant types compete with one another within but not among patches. Graminoids are known to have a significantly higher productivity (per unit area) than the other plant types others 2005, Johansson andothers 2006), whereas shrubs are typically least productive (Malmer and others 2005). The landscape-average NPP, provided as annual input to the model (see below), was thus partitioned among plant types reflecting their relative productivities, but maintaining the prescribed average NPP value (see equations (14 and B.1), respectively): where ANPP is the annual NPP at mean landscape level and ANPP x is the annual NPP for x plant type: shrubs, graminoids and three mosses, distributed according to D = 2.0, 1.5, 1.2, 1.0 and 1.0 for shrubs, graminoids, lawn, hummock and hollow mosses, respectively. ANPP is constant over time, because the ANPP x terms, affecting differential productivity through the D term in equation 14, are rescaled.

Permafrost/Freezing-Thawing Cycle
Freezing and thawing of peat soil is typical of arctic and subarctic conditions and leads to cryogenic landscape structures affecting plant productivity, decomposition and hydrological dynamics (Christensen and others 2004; Johansson and others 2006; Wania and others 2009b). To correctly calculate the fraction of ice and water in the peat soil, soil temperature at different depths must be estimated. We calculated the peat temperature at the centre of each layer (T L ) by adopting the simple analytical approach used in the LPJ-GUESS ecosystem model (Smith and others 2001;Sitch and others 2003): where T m is the mean of monthly mean temperatures for the last year (°C), EL is the exponential of oscillation lag, T s is the soil temperature at the centre of the peat layers, m and n are the regression parameters, AL is the oscillation lag in angular units at depth (from the surface to the centre of the peat layer) (see Table 1), and L is the conversion factor for oscillation lag from angular units to days (= 365/2p). Soil temperature is driven by surface air temperature which acts as the upper boundary condition.
The fraction of air, water and ice in each layer is updated daily based on the soil temperature in that layer, following the treatment of phase change described by (Wania and others 2009a). If the soil temperature calculated by equations 15-17 passes 0°C, ice content is set to zero and a corresponding amount of liquid water is added to the layer. Water and ice content are calculated by dividing the amount of frozen and melted water with total water-holding capacity. If a layer is totally frozen (100% ice), then it cannot hold additional water. In partially frozen soil, the sum of the fractions of water and ice is limited to water-holding capacity of that layer. The water and ice present in a particular peat layer influence soil thermal conductivity and heat capacity for that layer.

Study Area and Data Requirements
The model was applied based on the conditions at Stordalen, a subarctic mire in northern Sweden The warmest month is July, and the coldest month is February. The mean annual precipitation is low, but has recently increased from 304 mm  to 362 mm (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) (Johansson and others 2013). Overviews of the ecology and biogeochemistry of Stordalen are provided by Sonesson (1980), Malmer and others (2005) and Johansson and others (2006). Ecosystem respiration of Stordalen is lower relative to other northern peatlands due to low mean temperatures, a short frost-free season and the presence of discontinuous permafrost that keeps the thawed soil cooler and restricts decomposition rates (Lindroth and others 2007).
The model was forced with daily average air temperature and precipitation derived by combining millennium climate anomalies with the CRU climate dataset (Mitchell and Jones 2005). The method was explained in Chaudhary and others (2017a). From 1901-1912, the CRU TS 3.0 dataset of Mitchell and Jones (2005) and from the period 1913-2000 the observed dataset of Yang and others (2012) were used to force the model. The highspatial resolution (50 m), modern observed climate dataset was developed by Yang and others (2012) for the Stordalen site. In this dataset, the observations from the nearest weather stations and local observations were included to take into account the effects of the Torneträ sk Lake close to the Stordalen catchment on seasonal temperatures. The monthly precipitation data  for Stordalen at 50m resolution were downscaled from 10-min resolution using CRU TS 1.2 data (Mitchell and Jones 2005), a technique common for cold regions (Hanna and others 2005). The precipitation data were also corrected by including the influences of topography and by using historical measurements of precipitation from the Abisko Research Station record.
NPP To evaluate our model, simulation results were compared to an observed peat accumulation profile for the last 1000 years inferred from radioisotope dating of peat core sequences. The peat initiation started about 4700 years before present (BP) in the northern part and about 6000 BP in the southern part of the mire as a result of terrestrialization (Kokfelt and others 2010).

Simulation Protocol
To test hypotheses relating to peatland microtopography and peat accumulation, typical Stordalen initial conditions were maintained, where all the five plant types competed with each other based on the annual WTP with varying productivities and decomposition rates. The base model simulation (referred to as BAS in Table 3) was run for 1000 years with 50 patches. The model was initialized with a stochastically generated, uneven surface having sufficient water and ice, leading to a reasonably fixed WTP so that some areas remained drier, while some became completely saturated. Emulating the observed distribution of different microsites in the 1970s (Malmer and others 2005), the simulation was started from 1000 BP with almost 50% elevated patches, 30% patches representing intermediate areas and 20% representing depressions. This assumes no major shift in topographical, ecological and climate patterns prior to 1970.
A series of sensitivity tests was performed in order to determine the effects of key drivers and assumptions on the simulated peat accumulation and microtopographical dynamics (Table 3). In these experiments, the effects of ± 50% precipitation rates (P + 50 and P -50) and ± 5°C temperature (T + 5 and T -5) changes were simulated to ascertain the consequences of climate modifica-tions on peat accumulation and microtopographical changes. Implications of adjusting the homogenous species cover (HOM) and low surface roughness (LSR) on model outputs were also analysed. To isolate the effects of small-scale microtopography on peat accumulation, we removed the effects of source plant type on litter decomposition by running the model with only a single plant type (intermediate moss) in the HOM experiment. The LSR experiment enabled the influence of surface roughness on peat accumulation to be determined. Finally, to determine the effects of initial topographical structure on model output, in the RAN1 experiment the model was initialized repeatedly 50 times with randomly varying topographical structure.

Model Simulations
In the BAS experiment, the average total peat accumulation across 50 patches was 51.9 kg C m -2 in the last 1000 years BP, making a total cumulative peat depth increment of 49.5 cm (Figure 3 and Table 4). The simulated cumulative peat height is near to the observation-based estimate of 56 cm for the last 1000 years (Kokfelt and others 2010) and follows a similar trajectory. Areas dominated by hollow mosses (71.6 cm; 75.1 kg C m -2 y -1 ) and lawn mosses (64.9 cm; 65.5 kg C m -2 y -1 ) accumulated more peat relative to areas dominated by graminoids (58.2 cm; 61.1 kg C m -2 y -1 ), and dry mosses (31.1 cm; 32.6 kg C m -2 y -1 ) ( Figure 4A and B and Figure 5). The patches co-dominated by shrubs and hummocks mosses accumulated the least peat (25.6 cm; 26.9 kg C m -2 y -1 ). Although the initial decomposition rate of hollow mosses (0.1 y -1 ) and graminoids (0.25 y -1 ) was relatively higher, the saturated conditions in those patches limit daily decomposition rate leading to faster peat accumulation (Table 4). On the other hand, a fluctuating WTP increased the rate of decay in lawn mosses-covered ground despite a lower initial decomposition rate (0.05 y -1 ) and they accumulated less peat than patches dominated by hollow mosses. Shrubs and hummock mosses were assumed to have a moderate initial tissue decomposition rate (0.15 and 0.075 y -1 -see Table 1), yet showed the lowest peat accumulation in our simulations. This low rate of peat accumulation is the result of greater temperature-and dryness-driven rates of decomposition in elevated patches. Thus, the average rate of peat formation across the modelled landscape is ranked in the following order: Mho > Mim > Gr > Mhu > Sh ( Figure 5).
The mean annual simulated WTP fluctuated between -5 and -20 cm from the surface (Figure 4C). The mean annual simulated active layer depth (ALD) was shallow (20-30 cm) initially, but  reached 40-60 cm below the surface by the end of the simulation. The peatland's spatial heterogeneity, quantified by the standard deviation of patch height, declined over time ( Figure 4D). It is notable that graminoid-covered hollows deposited peat at faster rates compared to the elevated areas. Peat growth in intermediate areas was limited by high decomposition rates in elevated areas (a negative feedback). Therefore, we found that the spatial heterogeneity of the peatland decreased over the course of the simulation ( Figure 4D).
We found that many patches could be dominated by one or two plant types for a long period of time before reaching a threshold whereby an abrupt transition in dominance from one plant type to another occurred within a decade (not shown). On the other hand, in some patches, it was found that a new plant type could increase in abundance, after many years replacing the initial dominant plant types ( Figure 6A, B). The vegetation transition periods and FPC fluctuations differ between patches wherever this occurs, and during such tran-  Cyclicity between shrubs and hummock mosses was simulated in some elevated patches. The vegetation transition period lasts for 50-150 years during which time both types of vegetation can coexist ( Figure 6A). We found that a high rate of deposition of organic matter by hummock mosses could lead to a decline in moss fractional cover as the annual WTP drew down from the surface, resulting in a shift towards dominance by shrubs. However, the comparatively high decomposition rate of shrubs then slowed the growth of the peat column and the WTP again approached the surface in time, leading to suitable conditions for hummock mosses. A similar phenomenon results in a short cyclical transition between lawn and hummock mosses ( Figure 6B).

Sensitivity Analysis
Results from analysis of the model sensitivity to its forcing and initialization are shown in Table 5.
An increase in air temperature by 5°C (T + 5) resulted in thawing of permafrost and a lowering of the WTP owing to a higher evapotranspiration rate. Higher temperature also accelerated the microbial decomposition rate (see equation 5). This resulted  in a lower peat accumulation relative to the BAS experiment ( Figure 7). As decomposition increases due to lowering of WTP, a vegetation shift occurs and hummock mosses and shrubs completely occupy most of the patches and then peat accumulation declines, leading to WT approaching the surface once more. This enables the long-term resilience to the temperature change. A lowering of the surface air temperature by 5°C (T -5) substantially increased peat accumulation as a result of slower decomposition in cooler soils. The ALD became deeper, and patches hold less water (not shown). Evapotranspiration-driven water loss was decreased, further lowering the decomposition rate, but run-off rate increased as WTP neared the surface.
On decreasing the precipitation rate (P -50), we noticed a marginal decrease in peat accumulation because the peat column was almost frozen and only the upper 40-60 cm was active which can easily replenish even with such low levels of precipitation. Conversely, when the system is already saturated, any additional input of water (P + 50) will be removed because evaporation and surface run-off are increasing functions of WTP (see equations 8 and 9, respectively). However, a slightly higher average WTP in this simulation led to slightly higher peat deposition (Figure 7).
When heterogeneity between the patches was decreased (LSR) compared to the BAS experiment, the WTP approached the surface rather quickly leading to dominance by graminoids, hollow and lawn mosses and higher rates of peat accumulation (Figure 7). Over the course of the simulation, the patches dominated by hollow and lawn mosses increased more in height than graminoid-dominated areas, leading to a rougher surface by the end of the simulation.
In the HOM experiment, a homogenous cover of lawn moss accumulated 62.7 kg C m -2 (59.8 cm) of peat in 1000 years, which is high compared to the BAS experiment. When plant type-specific effects are removed, hydrological differences between microsites control the rate of peat accumulation and decomposition. Hollows added relatively more carbon than hummocks and intermediate areas due to their saturated conditions decreasing the overall decomposition rate. These differential rates of peat accumulation tended to cause convergence among patches, leading to a smoother peatland surface, in contrast to the LSR experiments described above. It can be seen from Figure A1 that the peat accumulation trajectory was quite sensitive to its initial microtopographical conditions (RAN1 experiments). This is because topography has a strong control on landscape WTP with cascading effects on plant distribution and peat accumulation and their distribution among patches within the landscape.

DISCUSSION
Classic peatland research focused mainly on the permanently saturated zone, the catotelm (Ingram  Figure 7. Peat accumulation over 1000 years in the sensitivity experiments (see Table 3 for details).
1982; Clymo 1984Clymo , 1992. The role of surface structures was highly generalized or ignored in earlier modelling studies. Our results demonstrate that surface structures may potentially play a critical role in determining the composition of the vegetation in different microsites depending upon the WTP, which in turn affects the plant litter quantity and quality, influencing the net rate of peat formation and overall peatland development.
Other recent studies have likewise emphasized the importance of surface microformations in peatland dynamics (Belyea and Malmer 2004;Belyea and Baird 2006;Pouliot and others 2011;Belyea 2013). The patterned surface creates a distinctive environment with unique plant cover, nutrient status, productivity and decomposition rates, constraining carbon fluxes at a range of scales, and affecting not just the peatland carbon balance but also its development.
Vegetation cover is important due to the consequences of plant litter quality on the accumulation rate of peat (Johnson and Damman 1991;Belyea 1996Belyea , 2013Thormann and Bayley 1997). Different plant species differ in productivity and produce litter that decomposes at different rates depending on its structural properties and chemical composition. This leads to highly variable spatial and temporal rates of carbon accumulation depending on vegetation type. Peat largely originating from graminoids and shrubs tends to decay faster than peat derived from Sphagnum mosses (Johnson and Damman 1991;Aerts and others 1999;Moore and others 2007;Strakova and others 2010). As a result, areas dominated by mosses accumulated comparatively more peat than the areas dominated by graminoids and shrubs in respective microforms (hummocks and hollows).
The BAS experiment highlighted the implications of plant litter quality on the net peat accumulation. Though graminoids are more productive than hollow mosses, their initial decomposition was also greater which led to relatively less peat accumulation in depressed sites. Likewise, areas codominated by shrubs and hummock mosses accumulated less peat than hummock moss areas despite shrubs being more productive (see Table 4). These results suggest that the inherent plant litter quality significantly affects the peatland growth and could be more important than the amount of litter deposited by the plant species.
Among mosses, Sphagnum species associated with hummock sites accumulated relatively less peat than hollow and lawn mosses due to a thicker aerated zone ( Figure 3A and 3B). Our simulation results suggest a hump-backed relationship be-tween the average rate of peat formation and water table position, and this is dictated by the differential WTP preferences assumed for the simulated plant types ( Figure 5). A hump-backed curve was similarly found to describe the relationship between the rate of peat formation and acrotelm depth at Ellergower Moss site by Belyea and Clymo (2001).
This hump-backed relationship was used to explain ecohydrological feedbacks and the nonlinearity that exists in the peatland system by (Belyea and Clymo 2001;Belyea 2013). There are two forces at work during the peatland microtopographical dynamics. One is stabilization force (negative feedback), and the other is a destabilization force (a positive feedback) (Belyea 2013), as shown in Figure 5. As a result of an external perturbation, a stable state is either pushed towards saturated (arrow I) or drier conditions (arrow II). The rising limb of the curve is considered unstable as any small climate perturbation may push the system towards drier and elevated path, but if the water saturation limit increases then the depressed patches may turn into deep-water pool, an irreversible stage (not represented in our model). On the other hand, a negative feedback (arrow IV) pushes back the system to the stable state by dampening the influence of external forcing and controlling the disproportionate growth of elevated microstructures. If this were not the case, the differences between the microforms would become exceptionally large. We allowed the peat to decompose in very dry conditions when the WTP drops below -40 cm and water content goes below 0.01 in the peat layers (see equation 4).
The total accumulated peat from the mixture of five plant types in the BAS experiment was lower than the HOM experiment, where only a single plant (lawn moss) cover was considered (Figure 7). The latter experiment was conducted to ascertain whether the decay rate is regulated by microhabitats or by plant species. We found that hollows accumulate more carbon when there are no differences in litter quality. Weltzin and others (2001) conducted an experiment in a homogeneously covered peatland and likewise observed higher peat accumulation in wet depressions. Our results suggest that the interaction of both initial decomposition rate and peat hydrology determine the amount of carbon deposited in patches.
Although empirical studies show that peatland microstructures may often remain stable for long periods of time, our findings demonstrate that although some areas (patches) may exhibit stability, overall a peatland can lose its heterogeneity in time ( Figure 4A and 4D). Furthermore, in some cases they may potentially also exhibit cyclic succession, as has often been discussed and hypothesized in early peatland vegetation studies (VonPost and Sernander 1910;Osvald 1923). In the BAS experiment, the cyclicity was simulated between the hummock mosses and shrubs, and between lawn and hummock mosses. They regenerated and replaced each other under a stable environment ( Figure 6A). In the former case (hummock mossesshrubs), the main reason for the cyclic transformation was the initial loss of surface roughness creating conditions suitable both for hummock mosses and for shrubs over higher areas. The elevated areas composed mainly of the peat derived from shrubs quickly lost peat mass as their decomposition rate was relatively high, forcing the peat column downwards (or to grow less quickly). This in turn created suitable conditions for hummock mosses as water neared the surface. When hummock mosses occupied these areas, they deposit enough carbon to cause the WTP to decline relative to the surface which again favoured dominance by shrubs. A similar cyclical phenomenon was noticed between lawn and hummock mosses ( Figure 6B), but soon disappeared as initial differences in patch height were reduced ( Figure 4D), resulting in higher WTP in all patches and a more favourable setting for lawn mosses. Tolonen (1985) showed continuous alternate streaks of dark and light layers from the peat sample of a raised bog in Maine and New Brunswick. These streaks were believed to correspond to lichen and Sphagnum peat, respectively, indicating that they are associated with hummocks and intermediate areas but not with hollows. Reflecting contrasting patterns of variability or stability among real peatlands, depending on hydromorphic factors, our results were generally consistent with the stable structure hypothesis (Tolonen 1971;Barber 1981;Seppa 2002), but also suggest that cyclic transformation may occur in intermediate and elevated areas under certain circumstances, or for transient periods. Such dynamics may be important for the evolution of large-scale peatland carbon balance under changing conditions such as climate warming, but may be challenging to predict with accuracy due to the strong influence of internal feedbacks suggested by our simulation results.
In general, heterogeneity was found to decrease over time in our simulations ( Figure 4D) and the peat landscape becomes progressively more dominated by Sphagnum species (Belyea and Malmer 2004). Hollow mosses dominate the deeper sites, whereas hummock mosses flourish in elevated areas. A stable state is attained where the Sphagnum species become increasingly dominant and occupy their own niche, as often seen in many peatlands (Belyea and Lancaster 2002), whereas the proportion of graminoid and shrub dominated patches decreases slowly over time ( Figure 4A). The decrease in aerated zone and dominance of Sphagnum species in turn promoted greater carbon accumulation and peat formation. We have also noticed that an abrupt transition from one dominant vegetation state to another state can lead to sudden shifts in the carbon accumulation trajectory (see grey lines in Figure 3).
In the BAS experiment, the vegetation composition in many peatland patches remained stable for a long period of time until the gradual climate forcing leads to a step-like shift in their vegetation states (not shown). Belyea and Malmer (2004) also noted that in many peatlands long periods of little or no change can end with a sudden, abrupt shift in vegetation composition and carbon balance. Subject to very strong forcing, these patches may follow a number of alternative pathways instead of following the dominant path, as can be seen (Figures 3 and 5) in our BAS and sensitivity experiments (T + 5 and T -5). Our results are in line with other studies (Hilbert and others 2000; Belyea 2013) and support the idea that successional dynamics can be nonlinear. The main identification of a nonlinear system is its sensitivity to initial conditions that in turn can cause it to follow number of alternative pathways, as reflected in our RAN1 experiment (Figure 7). Consistent with these findings, a regional compilation of palaeoecological records showed a diversity of histories of peatland development at different sites (Bunting and Warner 1998).
When the surface roughness was reduced (LSR), a total dominance of graminoids, hollow and lawn mosses was noticed as the annual WTP came closer to the surface on average. Over time, hollow and intermediate patches grew relative to graminoiddominated patches leading to a rougher peat surface, and by the end of the simulation some patches were occupied by hummock mosses indicating that system itself transforms, from smooth to heterogeneous surface and then stabilizes (Glaser and Janssens 1986;Pouliot and others 2011).
Peatland could lose or gain carbon in response to future climate warming. In our T + 5 experiment, higher temperatures lead to thawing of permafrost, higher soil temperatures and an increase of aerated (less waterlogged) soils due to higher evapotranspiration. These factors accelerate microbial decomposition (equation 5), leading to lower peat accumulation (Ise and others 2008; Dorrepaal and others 2009). As annual WTP drops, some patches come to be dominated by shrubs, which affect litter quality and peat accumulation. Studies have shown that NPP may also increase in some regions due to higher temperatures, a longer growing season and CO 2 fertilization (Euskirchen and others 2006;Chaudhary and others 2017b). In the absence of nutrient limitations, this could compensate for the loss of carbon from the system seen in our simulations, counteracting a potential positive feedback to climate warming (Beilman and others 2009;Jones and Yu 2010;Loisel and others 2012). Charman and others (2013) showed that variability in NPP has more influence than peat decomposition in determining peat accumulation patterns. Therefore, it is important to take into account both the dynamics of NPP in a warming climate and higher peat decay rates to determine future peat accumulation rates (Chaudhary and others 2017a). In this study we have focused on the interactions between vegetation dynamics and peat accumulation dynamics and have not taken into account the effect of climate change on the NPP to isolate these features and minimize the complexity of the model. Studies have also shown that colder and wetter conditions can decrease the peat decay rate due to an increase in anaerobic conditions and a lower substrate temperature, leading to higher peat deposition (Scanlon and Moore 2000;Frolking and others 2002;Malmer and Wallen 2004;Morris and Waddington 2011). This is also in line with our findings (T -5 simulation-see Figure 6).
In Stordalen, water table dynamics are largely controlled by underlying permafrost (Christensen and others 2004). In LPJ-GUESS, we have included soil freezing-thawing functionality which influences the model-generated data on NPP. The presence of frozen soil explains why no major change in annual WTP was noticed in the P ± 50 experiments. The majority of the lower peat layers were frozen, and the upper layers (surface 30-50 cm) were active which was quickly replenished under the restricted precipitation experiment (P -50). This in turn negligibly affected the peat accumulation relative to the BAS experiment. In the climate of Stordalen, soil freezing causes waterholding capacity of the patch to decrease with the result that the peatland requires less water to achieve a near-surface WTP compared to a nonpermafrost peatland. Conversely, in the P + 50 experiment, we found that adding extra water to peatlands which are already waterlogged with a high WTP (Figure 4) would not make much dif-ference to their WTP, vegetation cover and peat accumulation.
Some recent studies have pointed out that permafrost thawing may lead to an additional water input in low-lying and collapsed areas in high-latitude peatlands (Johansson and others 2006;Johansson and others 2013) under warming conditions. This may be expected to result in a shift in vegetation composition favouring graminoids and mosses, thereby affecting the carbon balance. In the Stordalen mire, permafrost underlying elevated areas is being degraded as a result of recent climate warming, with an increase in wet depressions modifying the overall carbon sink capacity of the mire (Christensen and others 2004;Malmer and others 2005;Johansson and others 2006). However, there is no soil subsidence functionality in our model, which explains in part why an increase in the wet patches was not simulated in the final years of the simulation.
Our results show that peatlands could exhibit different compositional and structural dynamics based on their initial topographical and climatic conditions and plant characteristics. This was also suggested by Strack and Waddington (2007), Ballantyne and others (2014), and Munir and others (2015) based on experimental treatments of peatland ecosystems, indicating a differential response of the observed microforms to both drought and warming, with the degree of the responses depending on the climatic region studied. These dynamics in turn will affect the net carbon balance from annual to millennial time scales. The microtopographical structure of peatlands may remain stable in some places, but overall its heterogeneity decreases if the surface is initially highly heterogeneous. Peatlands can also exhibit cyclicity in terms of vegetation and microtopographical structure under certain conditions (Figures 3 and 4). This implies that it may not be very easy to generalize the future behaviour of these microstructures in climate change impact studies and Earth system modelling (Baird and others 2013; Chaudhary and others 2017a, b). However, by highlighting the processes and parameters contributing the most variability in terms of carbon balance and vegetation cover, our parsimonious model approach represents a step towards the representation of such dynamics in regional vegetation and Earth system models, allowing their potentially important implications for regional and global carbon balance and biogeochemical feedbacks on the atmosphere to be explored and quantified.