Solid Earth Open Access Open Access Open Access

Abstract. Elevated nitrogen deposition and climate change alter the vegetation communities and carbon (C) and nitrogen (N) cycling in peatlands. To address this issue we developed a new process-oriented biogeochemical model (PEATBOG) for analyzing coupled carbon and nitrogen dynamics in northern peatlands. The model consists of four submodels, which simulate: (1) daily water table depth and depth profiles of soil moisture, temperature and oxygen levels; (2) competition among three plants functional types (PFTs), production and litter production of plants; (3) decomposition of peat; and (4) production, consumption, diffusion and export of dissolved C and N species in soil water. The model is novel in the integration of the C and N cycles, the explicit spatial resolution belowground, the consistent conceptualization of movement of water and solutes, the incorporation of stoichiometric controls on elemental fluxes and a consistent conceptualization of C and N reactivity in vegetation and soil organic matter. The model was evaluated for the Mer Bleue Bog, near Ottawa, Ontario, with regards to simulation of soil moisture and temperature and the most important processes in the C and N cycles. Model sensitivity was tested for nitrogen input, precipitation, and temperature, and the choices of the most uncertain parameters were justified. A simulation of nitrogen deposition over 40 yr demonstrates the advantages of the PEATBOG model in tracking biogeochemical effects and vegetation change in the ecosystem.


Introduction
Peatlands represent the largest terrestrial soil C pool and a significant N pool.Globally, peat stores 547 PgC (Yu et al., 2010) and 8 to 15 PgN, accounting for one-third of the terrestrial C and 9 % to 16 % of the soil organic N storage (Wieder and Vitt, 2006).Northern peatlands have accumulated 16 to 23 gC m −2 yr −1 throughout the Holocene and 0.42 gN m −2 yr −1 in the past 1000 yr on average (Vitt et al., 2000;Turunen et al., 2002;Limpens et al., 2006;van Bellen et al., 2011a, b).Carbon accumulation in peats has been primarily attributed to low decomposition rates, which compensate for the low production in comparison to other ecosystems (Coulson and Butterfield, 1978;Clymo, 1984).The two characteristic environmental conditions in northern peatlands' high water table (WT) and low temperature, play an essential role in preserving the large C pool by impeding material translocation and transformation in the permanently saturated zone (Clymo, 1984) Although the total N storage in peat is substantial, the scarcity of biologically available N induces a conservative manner of N cycling in peatlands (Rosswall and Granhall, 1980;Urban et al., 1988).Sphagnum mosses are highly adapted to the nutrient-poor environment and successfully compete with vascular plants through a series of competition strategies, such as inception of N that is deposited from the atmosphere, internal recycling of N, and a minimized N release from litter with low decomposability (Damman, 1988;Aldous, 2002).
Climate change and elevated N deposition are likely to alter the structure and functioning of peatlands through interactive ways that are incompletely understood.In general, drought and a warmer environment were found to affect Published by Copernicus Publications on behalf of the European Geosciences Union.
Y. Wu and C. Blodau: PEATBOG vegetation composition by suppressing Sphagnum mosses and promoting vascular plants (Weltzin et al., 2003), which in turn alters litter quality, C and N mineralization rates (Keller et al., 2004;Bayley et al., 2005;Breeuwer et al., 2008), and the C and N balance (Moore et al., 1998;Malmer et al., 2005).In northern peatlands, nitrogen is often a limiting nutrient and regulates the rates of C and N cycling and individual processes, and thus also controls elemental effluxes to the atmosphere and discharging streams.Excessive N entering peatlands could induce changes in various processes that may lead to non-linear and even contrasting consequences with respect to C and N budgets, especially on longer timescales.For example, experimentally added N was found to increase photosynthetic capacity and growth of several Sphagnum species up to ca. 1.5 gN m −2 yr −1 before causing their decline at low N background sites (Williams and Silcock, 1997;Granath et al., 2009).However, at high N background sites such effects occurred up to 4 gN m −2 yr −1 (Limpens and Berendse, 2003), which raises the question of how peatland ecosystems adjust their structure and functioning to long-term N deposition.Survey studies across N deposition gradients ranging from 0.2 to 2 gN m −2 yr −1 demonstrated a relation between N deposition and litter decomposition rates (Bragazza et al., 2006), in addition the effects seemed to depend on litter quality (Bragazza et al., 2009;Currey et al., 2010) and deposited N forms (Currey et al., 2010).In both long-term N fertilization experiments and survey studies, an increase in N content in the surface peat and in the soil water was observed at the high N sites (Xing et al., 2010) but enhanced N effluxes in the form of N 2 O remained elusive (Bubier et al., 2007).In contrast, N 2 O emission was found in short-term N and P fertilization experiments (Lund et al., 2009).Laboratory and field experiments aiming to quantify the combined effects of temperature, WT and N elevation have thus often arrived at contradictory conclusions, due to the interplay of effects in time and space (Norby et al., 2001;Breeuwer et al., 2008;Robroek et al., 2009).Furthermore, elevated N deposition was recently suggested to affect soil temperature and moisture through changes in the vegetation community with potential feedbacks on elemental cycles (Wendel et al., 2011).
Ecosystem modeling has become an important approach in analyzing the interacting effects of climate and N deposition on peatlands and in making long-term predictions; examples are provided by PCARS (Frolking et al., 2002), ecosys (Dimitrov et al., 2011), Wetland-DNDC (Zhang, 2002), and MWM (St-Hilaire et al., 2010).While models have been thoroughly developed to investigate peatland C cycling (e.g.PCARS, MWM), there have been few attempts to integrate N cycling into peatland models, although N is mostly considered to be the limiting factor on primary production (Heijmans et al., 2008).In the mentioned models, N is generally passively bound to C pools by C/N ratios, while active nitrogen transformation and translocation among N pools is omitted.
To make progress towards closing this gap, we present a novel model for the analysis of the coupled C and N cycles in northern peatlands.The model is designed to fulfill the following objectives: (1) to clarify the interaction between C and N cycling in vegetation, soil organic matter and soil water; (2) to determine key processes that control the C and N balance of northern peatlands in the short and long-term; (3) to quantify C and N pools and cycling rates in peatlands; (4) to characterize their sensitivity to N availability and climate change; and (5) to predict the combined impact of elevated N deposition and climate change on peatland C and N cycling.
In this paper, we focus on the integration of C and N cycling through vegetation, soil organic matter and soil water, the coupling of C and N throughout the ecosystem, and the consistency of mass movements between pools.We first highlight the structural design and principles that governed the modeling process, and then explain the components of the model by focusing on the individual submodels.To improve readability of the text the equations are listed in the appendix.We subsequently present an evaluation of the simulated WT dynamics, C fluxes, depth profiles of CO 2 and CH 4 in soil water, and C and N budgets.The model output is compared against observations for the well-characterized Mer Bleue Bog (MB), Ontario, Canada.We also present sensitivity analyses for environmental controls, such as temperature, precipitation, and N deposition, and for some calibrated key parameters.Finally we demonstrate the potential of the model for analyzing the effects of experimental long-term N deposition and climate change.

Model description
The PEATBOG (Pollution, Precipitation and Temperature impacts on peatland Biodiversity and Biogeochemistry; see acknowledgements) model version 1.0 was implemented in Stella®and integrates four submodels: environment, vegetation, soil organic matter (SOM), dissolved C and N (Fig. 1).The environment submodel generates daily WT depth from a modified mixed mire water and heat (MMWH) model (Granberg et al., 1999) and depth profiles of soil moisture, peat temperature and oxygen concentration.The vegetation submodel simulates the C and N flows and the competition for light and nutrients among three plant functional types (PFTs): mosses, graminoids and shrubs.Most of the algorithms of plant physiology were adopted from the Hurley pasture (HPM) model (Thornley and Verberne, 1989;Thornley et al., 1995;Thornley, 1998a).Modifications were made for mosses and for the competition among PFTs in the nutrient-poor environment.Litter and exudates from the vegetation submodel flow into the SOM submodel and are decomposed into dissolved C and N. The dissolved C and N submodel tracks the fate of dissolved C and N as DOC, CH 4 , CO 2 and DON, NH + 4 , and NO − 3 .The model does not consider hummock-hollow microtopography of peatlands, which in other studies had no statistically significant effect when simulating ecosystem level CO 2 exchange (Wu et al., 2011).

Model structure and principles
The following three principles were imbedded in the model in terms of scale, resolution and structure:

High spatial and moderate temporal resolution
In comparison to other biogeochemical process models of peatland C cycling (Frolking et al., 2002;St-Hilaire et al., 2010) that primarily focus on the ecosystem-atmosphere interactions, we increased the vertical spatial representation and kept the temporal resolution fairly low.We divided the belowground peat into 20 layers (i) with a vertical resolution of 5 cm except for an unconfined bottom layer.This structure applies to all belowground pools and processes.The rationale for the comparatively fine spatial resolution lies in the critical role of soil hydrology for the C and N cycles and the necessity to represent physical and microbial processes (Trumbore and Harden, 1997).Spatial distributions of water and dissolved chemical species are generated and mass movement and balances are examined throughout layers and pools, which allows for tracing the fate of C and N belowground.The high resolution allows to explicitly include the activity of plant roots and their local impact on C and N pools.Plant roots showed morphological changes upon WT fluctuation and nutrient input in bogs (Murphy et al., 2009;Murphy and Moore, 2010).Root litter also provides highly decomposable organic matter to deeper peat and serves as a substrate for microbial respiration.Moreover, roots can act as sensitive conductors of N deposition to deep peat via root chemistry and litter quality (Bubier et al., 2011;Bragazza et al., 2012).The layered structure assists in mapping the belowground micro-environment for simulating the sensitive interactions of soil moisture, roots and microbial activity.The model computes and simulates processes on a daily time step, as does for example the HPM model (Thornley et al., 1995) and the wetland-DNDC model (Zhang, 2002).The moderate temporal resolution is adequate for the model soil C in the short-and long term (Trettin et al., 2001).

Stoichiometry controls C and N cycles
We did not stipulate critical mass fluxes as constraints on C and N cycling.Instead these constraints are generated in the model from changes in biological stoichiometry.This structure has the advantage that the interactions between C and N fluxes and temporal and spatial changes in pools sizes control the mobility of the elements.As in some terrestrial C and N models (Zhang et al., 2005), N flows are driven by C/N ratio gradients from low C/N ratio to high C/N ratio compartments.The C/N ratios of all pools are in turn modified by their associated flows, reflecting the organisms' requirement to maintain their chemical composition in certain ranges.Results from field manipulation experiments suggested thresholds of the N deposition level, above which the Sphagnum moss filter fails and mineral N enters soil water (Lamers et al., 2001;Bragazza et al., 2004).Flux-based critical loads of N for Sphagnum moss were suggested as the high end of the Sphagnum tolerance range, where the values are between 0.6 gN m −2 yr −1 (Nordin et al., 2005) and 1.5 gN m −2 yr −1 (Vitt et al., 2003).Threshold values in stoichiometry terms appear to be less variable, ranging from 15 mgN g −1 (Van Der Heijden et al., 2001;Xing et al., 2010) to 20 mgN g −1 dry mass (Berendse et al., 2001;Granath et al., 2009).The critical load of ca. 1 gN m −2 yr −1 was linked to a stoichiometry thresholds of 30 (N/P ratio) and 3 (N/K ratio) in Sphagnum mosses (Bragazza et al., 2004).The model internally generates C/N ratios, or C/N/P ratios, for all compartments to control the N flows in plants and microorganisms.

Consistent conceptualization of carbon and nitrogen reactivity
Differences in the mobility of C and N compartments were implemented using a two-pool concept throughout the model.Similar to decomposition models that distinguish the quality of soil organic matter (Grant et al., 1993;Parton et al., 1993), C and N are presented in labile (L) and recalcitrant (R) pools in SOM.In addition, the model differentiated C and N pools based on quality in vegetation, into structural (struc) pools (Fig. 2).The pasture vegetation model HPM (Thornley et al., 1995;Thornley, 1998b) was adopted, where C and N in grass and legumes were separated in structural and substrate pools in shoots (sh) and roots (rt) for 4 age categories.Considering our focus on competition between plant functional types, vegetation was not conceptualized in terms of age categories but instead classified into 3 plant functional types (PFTs) (j : 1 = mosses, 2 = graminoids and 3 = shrubs) that are characterized by distinctive ecological functions (Fig. 3) in our model.The plant functional types differ in the decomposability of the litter, which was represented by the different mass fractions of the labile carbon pool in the litter.The fraction of labile litter was assumed to be 0.1, 0.3 and 0.2 in mosses, graminoids and shrubs, respectively (Inglett et al., 2012).Once the litter is deposited the litter merges into one labile and one recalcitrant soil organic matter pool.The remaining fraction of the plant litter is assigned to be recalcitrant and represents the input into the recalcitrant soil organic matter pools.Thus, the composition of plants, as a result of net primary production and litter fall, is adjusted to physical conditions and N input and alters SOM quality via changes in litter quality (Q).

Structural adaptations for modeling peatland biogeochemistry
Modifications were made to the adopted algorithms of the MMWH and HPM models for compatibility with our modeling purpose and model structure.The main modifications and novel features of the PEATBOG model are:

Competition among Plant Functional Types (PFTs)
Plant functional types compete for light and nutrients through their morphology and nutrient utilization.We modified the algorithms of competition among plant functional types for these controls to better represent the shading effects among PFTs and the nutrient-poor environment.Competition among plants was modeled using PFTs previously, where the depth and biomass of roots mainly determined superiority in competition (Van Oene et al., 1999;Pastor et al., 2002;Heijmans et al., 2008).We focused instead on the effect of light for PFT competition that is controlled by shading effects through canopy layers (Fig. 3).This differs from the utilization of the leaf area index, which determines the share of total photosynthesis in the HPM model (Thornley et al., 1995).In the PEATBOG model, the uptake of N is also modified to be specific for each soil layer and PFT.It includes the uptake of three forms of N in the PFTs so that N availability varies for roots of each PFT in the same location.In addition to inorganic N sources (NH + 4 and NO − 3 ), as modeled in some C and N cycling models (Aber et al., 1997;Van der Peijl and Verhoeven, 1999), DON is included as a third N source, acknowledging its abundance (Moore et al., 2005) and potential importance in nutrient-poor environments, such as bogs (Jones et al., 2005;Nasholm et al., 2009) (Fig. 3).

Decoupling of O 2 boundary and WT boundary
The interface between oxic and anoxic conditions and unsaturated and saturated peat (i.e. the water table position, WT) are separately modeled and control biogeochemical and physical processes, respectively.Recent findings questioned that the long-term WT is the sole control on biogeochemical processes in peat as well as the acrotelm and catotelm concept in modeling of peatlands (Morris et al., 2011).Meanwhile O 2 was found well above and below the WT in peats, for instance during drying and rewetting experiments in a degraded fen site with dense soil (Estop-Aragonés et al., 2012).The decoupling of redox conditions from the WT spatially and temporally in dense soils is potentially important for the partitioning of respired C into CO 2 and CH 4 during the decomposition of peat.We calculated O 2 concentration in each layer to regulate energy-limited processes such as CH 4 oxidation and peat decomposition.Water table, on the other hand, serves as a control on moisture-limited biological or physical processes, such as root metabolism and diffusion.The belowground controls on CH 4 production and emissions and the advantages and disadvantages of our representation of oxygen and soil moisture dynamics will be further discussed in a future manuscript.

Change in Figures
Please change Fig. 2 to the following:

Submodel 1 -environmental controls
Physical boundary conditions, such as day length, degree days, water table depth, soil moisture, temperature and depth profiles of O 2 , are generated by the model to control physiochemical and biological processes.Day length (DL), which in the model controls photosynthesis, varies for geographic position of the site and day of year.The daily day length value is obtained from the angle between the setting sun and the south point, which in turn is calculated from the declination of the earth and the geographical position of the site (Brock, 1981) (Appendix A, Eqs.A14 and A15).Declination of the earth is the angular distance at solar noon between the sun and the equator and positive for the Northern Hemisphere.The value of declination is approximately calculated by Cooper (1969) using the day of the year.
Temperature is modeled by sinusoidal equations (Carslaw and Jaeger, 1959) and modified by converting a dampening depth into thermal conductivity (Appendix A, Eq.A13).Thermal conductivity (K thermal ) is adjusted for each layer for peat compaction and snow coverage that delays the thermal exchange in winter and early spring (Fig. S1a, Supplement).
Degree days (DD) represent the accumulation of cold days and trigger defoliation (Frolking et al., 2002;Zhang, 2002).Similar to other models, defoliation occurs on the day when DD reaches minus 25 degrees, with accumulated temperature of lower than 0 degrees after day 181 of the year (1 July in non-leap years).
Water table (WT) depth is simulated by calculating the water table depth from the water storage of peat using a modified version of the Mixed Water and Heat model (MMWH) (Granberg et al., 1999).Precipitation and snowmelt represent water inputs, and are obtained from local meteorological records, instead of modeling the snow cover.Evapotranspiration (EPT) is the water output from the peat and vegetation surface via evaporation and transpiration, which are regulated by temperature and vegetation characteristics.Different from the authors' original approach, the EPT rate per unit of the peatland surface is calculated from a base EPT rate and multipliers of plant leaf area (Reimer, 2001) (Appendix A, Eq.A3), daily air temperature (Fig. S1b), daily average photosynthetic active radiation (PAR), and a factor of WTD and rooting depth (Lafleur et al., 2005a) (Fig. S2c).A maximum water storage was added to allow overflow once the WT rises above the peat surface.WTD is then obtained from linear functions of water storage as in the MMWH model but with depth-dependent slopes (Appendix A, Eq.A8).The WT layer is defined as the layer in which the WT is located.
Depth profiles of soil moisture (m 3 water m −3 pore space) are generated by the Van Genuchten's soil water retention equation, parameterized by Letts et al. (2000) for peatlands (Appendix A, Eq.A9).Porosity is a function of depth derived from field measurements for the Mer Bleue Bog (Blodau and Moore, 2002).
In order to simulate exports of dissolved C and N without modeling water movement explicitly, runoff was distributed over 20 layers and divided into horizontal and vertical flows (Fig. 4, Appendix A, Eq.A4-A7).The vertical advection rate depends on slope and is determined as a fraction of the total runoff.It is consistently applied to all layers.The remaining runoff is horizontally distributed among layers according to the vertical hydraulic conductivity distribution.In the Mer Bleue Bog, saturated hydraulic conductivity rapidly declines with depth in the acrotelm, ranging from 10 −7 to 10 −3 m s −1 and reaching 10 −8 to 10 −6 m s −1 in the catotelm (Fraser et al., 2001b).In layers above the WT, the actual hydraulic   conductivity is lower when pores are unsaturated (Hemond and Fechner-Levy, 2000) (Fig. S1d).
The depth profiles of O 2 concentrations are simulated to locate the oxic-anoxic interface.Oxygen diffuses from the surface to deeper soil layers and is consumed directly or indirectly by the oxidization of peat C to CO 2 (Appendix A, Eq.A12).For the simulation of oxygen-dependent biogeochemical processes we chose a dichotomous distribution of O 2 , where the boundary of oxic/anoxic conditions is set at 5 µmol L −1 (Liou et al., 2008).

Submodel 2 -vegetation
Carbon in vascular plants is represented by four pools: shoot substrate C (sh subsC), root substrate C (rt subsC), shoot structural C (sh strucC), and root substrate C (rt subsC) (Fig. 2).Substrate C and structural C refer to metabolic activated C and recalcitrant C, respectively.Substrate pools conduct metabolic activities (i.e.photosynthesis, respiration) and structural pools perform phenological activities (i.e.growth, litter production).The flow from substrate C to structural C leads to plant growth (Appendix A, Eq.A33).Each C pool or flow is bound to an N pool or flow by the C/N ratio of the specific pool.Furthermore, shoots are divided into stems and leaves and roots into coarse and fine roots by ratios specific to the PFT.Mosses are represented by 4 aboveground pools and two compartments: capitulum and stem.The C and N contained in exudates are transferred from the vegetation into the uppermost labile C and N pools in the soil.Unlike N uptake by vascular plants from soil water, N uptake by mosses is restricted to atmospheric supply.
Most C and N material flows are driven by C concentration gradients except for a few processes controlled by N (i.e.N uptake, N recycling from litter production).The phenology and competing strategies of PFTs are modeled as follows: (1) considering the seasonal C and N loss in leaves of deciduous shrubs; (2) PFT-specific N flows during growth, recycling and litter production; (3) competition among PFT is implemented through shading effects, tolerance to moisture and temperature, distribution of C and N among shoots and roots, as well as turnover rates.In general, the photosynthetic nutrient-use efficiency (the ratio of photosynthesis rate and nitrogen content per leaf area) is higher in herbaceous than in evergreen woody species (Hikosaka, 2004).The growth rates in deciduous species (graminoids and deciduous shrubs) are higher than in evergreen shrubs, which in turn is higher than in mosses (Chapin III and Shaver, 1989).Graminoids are more competitive in the deep soil attributed to the longer roots (Murphy et al., 2009).Mosses have the advantage of aboveground N uptake and filtration.Below we discuss the modeling of these competition strategies.

Photosynthesis (PSN) and competition for light
Competition for PAR is implemented through shading effects.The light level that reaches a specific PFT after interception by a taller PFT determines the C assimilation of this PFT (Fig. 3).For each PFT, canopy PSN is integrated from daily leaf PSN by a light attenuation coefficient (k ext ), leaf area index (LAI) and day length (DL) (Appendix A, Eq.A38).The coefficient k ext is unitless, the values are 0.5 for graminoids (Heijmans et al., 2008), 0.97 for shrubs (Aubin et al., 2000), and assumed to be 0.9 for mosses.LAI is determined by leaf structural C mass and specific leaf area (SLA) of the PFT.The PSN rate for the top canopy layer of each PFT (LeafPSN j ) is calculated by a non-rectangular hyperbola (Fig. S2f, Appendix A, Eq.A40).The two parameters α j and ξ control the shape of the hyperbola curves.Parameter α j represents the photosynthetic efficiency, which is controlled by WT depth, the air temperature (T air ) and atmospheric CO 2 level (CO 2,air ) (Appendix A, Eq.A42).The spring PSN of mosses starts when the snow depth falls below 0.2 cm.The variable LI j is the PAR incepted by the canopy of PFT j (umol m −2 s −1 ).The assumptions here were that radiation diminishes along with canopy depth and each canopy depth contains one PFT solely.
The asymptote of leaf photosynthesis rate (P max in gCO 2 m −2 s −1 ) is regulated by T air , CO 2,air , WT depth, N content in plant shoots and the season.The maximum PSN rate (P max,20 , g CO 2 m −2 s −1 ) occurs in an optimal environment, is also referred to as PSN capacity, and is often derived from measurements.The values of P max,20 vary among and within growth forms and follow the general sequence of deciduous > evergreens > mosses (Chapin III and Shaver, 1989;Ellsworth et al., 2004).The maximum PSN rate P max,20 is 0.002 g CO 2 m −2 s −1 for graminoids and mosses following HPM (Thornley, 1998a), and 0.005 g CO 2 m −2 s −1 for shrubs based on the ranges in Small (1972).The temperature dependence (f T ,P max,j ) of P max is conceptualized as sigmoidal curve with PFT-specific optimal, maximum and minimum temperature for photosynthesis and curvature q (Fig. S2e, Appendix A, Eq.A43).The WT depth dependency of P max (f m,P max,j ) for mosses follows Frolking et al. (2002) and is an exponential function with PFT-specific base (a w,j ) for vascular plants (Fig. S2a and b).The model considers season and nutrient availability effects on P max .Seasonal change (f season,P max ) affects mosses alone between 0 to 1 and was derived from the maximum rates of carboxylation (V max ) in spring summer and autumn (Williams and Flanagan, 1998) (Fig. S2c).
Potential N stress on photosynthesis is modeled by using PFT-specific photosynthetic N use efficiencies.Although there are interacting controls on the N economy of plant photosynthesis, such as N effects on Rubisco activity, Rubisco regeneration and the distribution of N in leaves, there seems to be a generalized linear relation of foliar N content and PSN capacity across growth forms and seasons (Sage and Pearcy, 1987;Reich et al., 1995;Yasumura et al., 2006).The ratio of PSN capacity and foliar N concentration is defined as photosynthetic nitrogen use efficiency (PNUE) (Field and Mooney, 1986).In general, evergreens have lower PNUE and larger interception than the deciduous shrubs (Fig. S2d, Appendix A, Eq.A47) (Hikosaka, 2004).To reflect N use strategies of growth forms, we implemented PNUE values for PFTs following the sequence: graminoids > shrubs > mosses, and interception values reversely.In addition, a toxic effect (f N,toxic ) is applied with regard to mosses when the substrate N concentration exceeds the maximum N concentration at 20 mg g −1 (Granath et al., 2009).

Competition for nutrients
PFTs compete for N through two processes: filtration of deposited N by mosses and the uptake of N among vascular plants roots.Nitrogen deposited from the atmosphere is first absorbed by moss and then enters soil water to become available to vascular plants.The N/P ratio of mosses is used as a regulator of N pathways and an indicator of N saturation in mosses.A fraction of 95 % of the deposited N is absorbed by moss until the N/P ratio reaches 15 (Aerts et al., 1992), above which N absorption decreases owing to the co-limitation of N and P on PSN rates.We assume mosses become N saturated when the N/P ratio exceeds 30 (Bragazza et al., 2004), above which the uptake fraction declines to zero.Due to the lack of P pools in the current model version, the initial moss N/P ratio is assumed to be 10 in mosses (Jauhiainen et al., 1998).
Graminoids have a larger rt k than shrubs, indicating more roots in deeper layers that allow utilization of N in deeper peat.The N uptake rate is affected by the surface area rather than the biomass of the fine roots.Specific root lengths LV j that vary with root diameters are used to convert the dry biomass to the surface area of roots (Kirk and Kronzucker, 2005).The diameters of the fine roots were set to be between 0.005 to 0.1 cm for the "true fine roots" that are responsible for N uptake (Valenzuela-Estrada et al., 2008).Nitrogen uptake is modeled using Michaelis-Menten equations (Appendix A, Eq.A71-A73), controlled by the soil temperature, the root biomass of the layer and the substrate C and N concentrations in plants.Parameters V max and K m for the DIN uptake were derived from the model of Kirk and Kronzucker (2005) while those for DON uptake were calibrated based on one of the few quantitative studies for an Arctic Tundra (Kielland, 1994), where V max for DON uptake was 0.0288 to 0.048 mmol g −1 day −1 for shrubs (Ledum) and 0.012 to 0.096 mmol g −1 day −1 for graminoids (Carex/Eriophorum).The effects of substrate C and N concentration in plants on N uptake rates were derived from the HPM model (Thornley and Cannell, 1992).The half saturation constant of substrate N was adjusted to be smaller for shrubs and mosses than for graminoids.The temperature influence on N uptake is modeled using Q 10 functions for active NO − 3 uptake and linear functions for passive NH + 4 uptake (Glass et al., 2001;Williams and Miller, 2001;Miller and Cramer, 2004).Despite the abundance of DON in soil water, which is one order of magnitude larger than the concentration of DIN in the field (Kranabetter et al., 2007;Nasholm et al., 2009), the capability of DON uptake by plants is limited to low molecular weight DON (e.g.glycine, aspartate and glutamate) (Jones et al., 2005).We assumed a fraction of 0.2 of total DON concentration to be bio-available to plants, according to reports on arctic tundra and two permafrost taiga forests (Jones and Kielland, 2002;Atkin, 2006).Pools of NH + 4 , NO − 3 , and DON are simulated in the dissolved C and N submodel.

Submodel 3 -soil organic matter dynamics
The soil organic matter (SOM) submodel simulates peat decomposition and accumulation using a multi-layer approach.The litter produced from the vegetation submodel is added to the topsoil layer and into the rooted layers of the peat.In each layer, C and N are present in labile (L) and recalcitrant (R) pools.The decomposition of each SOM pool was modeled following the single pool model of Manzoni et al. (2010).Pool L and R are decomposed simultaneously at rates that are determined by their C/N ratios, an environmentally controlled decomposition rate constant k, and the availability of mineral N. Three fates of the decomposition products are possible: (1) leaching as dissolved organic matter (DOM), (2) re-immobilization into microbial biomass, and (3) conversion into dissolved inorganic carbon (DIC) and dissolved inorganic nitrogen (DIN).DOM was extracted from SOM pools by a constant fraction, which is empirically related to the local precipitation level of the site (Appendix A, Eqs.A90 and A96).The value used here (0.05) is slightly smaller than the lower end (0.06) of the suggested range for ecosystems in general (Manzoni et al., 2010), owing to the small hydraulic conductivity in northern peatlands.The remaining SOM is either mineralized into dissolved inorganic matter or immobilized into microbial biomass with a microbial efficiency (e), indicating the immobilized fraction of the decomposed SOM (Appendix A, Eq.A84).Parameter e is empirically calculated from the initial C/N ratios of the SOM pools, which in turn is controlled by the composition of litter produced from each PFT.For simplicity, microbial biomass is considered as a constant part of SOM.The actual N decomposition rate, excluding for the N immobilization to microbial biomass, can be either positive or negative.Positive rates reveal net mineralization from SOM N pools to dissolved NH + 4 pools and negative rates indicate net immobilization.The "critical N level" is used as an indicator of the N concentration at which immobilization balances mineralization (Berg and Staaf, 1981).The "critical N level" varies according to the C/N ratio of microorganisms, the DOM leaching fraction, e and another factor representing the N preferences of microorganisms during decomposition (α E N prefer ) (Appendix A, Eq.A86).The nitrogen preference of microorganisms (α E N prefer ) is a multiplier larger than 1 and is limited by the asymptotic C/N ratio of SOM at decomposition equilibrium (Appendix A, Eq.A95).
In addition to the control of N concentration in SOM, the availability of soil mineral N also affects the decomposition rates.Nitrogen addition experiments showed neutral or negative effects on the decomposition rates of SOM due to contrary effects on the decomposition of labile and recalcitrant OM: a decrease in the decomposition rates of more recalcitrant OM and an increase in that of more labile OM (Neff et al., 2002;Janssens et al., 2010;Currey et al., 2011).We adopted the quantitative relation from the Integrated Biosphere Simulator model (IBIS) (Liu et al., 2005), by converting mineral N contents to DIN concentrations in each layer (Fig. S3d).Nitrogen mineralization is inhibited while N immobilization is promoted by increasing DIN concentration up to 200 µmol L −1 .The decomposition rate constants k are regulated by substrate quality (q), soil moisture (fm dec ), soil temperature (f T dec ) and inhibition factors accounting for the decrease in Gibbs free energy due to the accumulation of end products (i.e.CO 2 , CH 4 ) in the saturated soils (Appendix A, Eq.A87).The decrease in k with depth is modeled based on the "peat inactivation concept" (Blodau et al., 2011) rather than only linked to anoxia (Frolking et al., 2002) or redox potential (Zhang, 2002), as in other models.The essential idea of this concept is that the transport rate of decomposition products controls the decomposition rate in the saturated anoxic soils (Fig. S3).The inhibitions factors are values between 0 and 1 based on CO 2 and CH 4 concentrations according to the inverse modeling results in Blodau et al. (2011) (Fig. S3a and b).
The intrinsic decomposability of the substrate (L or R) determines the base decomposition rate constant (k Cpot ).Due to the conceptual inconsistency of k Cpot in experiments (Updegraff et al., 1995;Bridgham et al., 1998), we calibrated the values of k Cpot from the long-term simulations in the spin-up runs.The moisture and temperature effect on the decomposition is each pool is modeled similar to the PCARS model (Frolking et al., 2002), with the Q 10 value of the decomposition of L pools (2.3) smaller than of that of R pools (3.3) (Conant et al., 2008(Conant et al., , 2010)).

Submodel 4 -dissolved C and N
The model contains 3 dissolved C pools: CH 4 , CO 2 and DOC and 4 dissolved N pools: NH + 4 , NO − 3 , NO − 2 and DON in each belowground layer (Fig. 2).Because decomposition proceeds and is controlled through the SOM pools, DOC and DON are considered to be an end product, and are only removed by runoff.The production of DOC, DIC, DON and NH + 4 are inputs from the SOM and the vegetation submodels.The production of DIC is further partitioned into the production of CH 4 and CO 2 in the anoxic layers.
The partitioning of respired C into CO 2 and CH 4 in the saturated layers depends on the presence of alternative electron acceptors (i.e.SO 2− 4 , NO − 3 and likely humic substances) for the terminal electron accepting processes (TEAP) (Conrad, 1999;Lovley and Coates, 2000).In previous studies, the ratio of CO 2 /CH 4 production and the production rates of CH 4 was modeled as a function of WT depth (Potter, 1997;Zhuang, 2004), or by microbial activities using Michaelis-Menten kinetics (Segers and Kengen, 1998;Lopes et al., 2011).Following the concept put forward by Blodau (2011), we modeled the CH 4 production rate by an energy-limited Michaelis-Menten kinetics.
We built an equation group based on the valance balance of the overall oxidation-reduction process and the mass balance of C (Appendix A, Eq.A121).The first equation (Appendix A, Eq.A121) denotes that CO 2 and CH 4 are the only inorganic C products (DIC) from the decomposition of SOM.The second equation was deduced from the valance balance of CO 2 (+4) production and CH 4 (−4) production from organic C, assuming an initial oxidation state of zero as found in carbohydrates.The production of CO 2 (CO 2 pro i ) is the result of the stoichiometric release of CH 4 (CH 4 pro i ) from fermentation and subsequent methanogenesis, and the consumption of electron acceptors (CO 2 pro EA,i ) in units of electron equivalents.The acronym EA represents electron acceptors other than CO 2 , including NO − 3 , SO 2− 4 , and humic substances (HS).
In anaerobic systems, electron acceptors are consumed by terminal electron accepting processes that competitively consume H 2 or acetate.Individual processes predominate according to their respective Gibbs free energy gain, usually in the sequence NO − 3 , Fe (III), humic substances (HS), SO 2− 4 and CO 2 (Conrad, 1999;Blodau, 2011).Owing to the extremely fast turnover of H 2 pools in peat, the Michaelis-Menten approach is not suitable for modeling CH 4 production in models running on a daily time step when H 2 is considered the substrate.To avoid modeling the pools of H 2 and acetate explicitly, the current model with daily time step focuses on the electron flow from complex organic matter to all TEAPs, instead of modeling each microbial process explicitly.In ombrotrophic systems like bogs, only SO 2− 4 , NO − 3 and HS are considered relevant electron acceptors.The CO 2 production from SO 2− 4 and NO − 3 reduction are calculated from the valance relations (Appendix A, Eq.A122), One mole of SO 2− 4 being reduced to HS − provides 8 mole of electrons (S(+6) → S(−2)) and 1 mole of NO − 3 release 5, 4 and 3 moles of electrons when being reduced to NO, Humic substances have recently also been identified as electron acceptors (Lovley et al., 1996;Heitmann et al., 2007;Keller et al., 2009) and require some consideration.Reduction of humic substances may be a significant CO 2 source in anoxic peat, where a large fraction of the total CO 2 production typically cannot be explained by consumption of known electron acceptors (Vile et al., 2003b).Although peat stores a large amount of organic carbon as humics, likely only a small fraction of it is redox active (Roden et al., 2010).The redox-active moieties in humics have been identified as quinones, here called DOM-Q (Scott et al., 1998).Electron accepting rate constants of HS in sediments were reported to be 0.34 h −1 and 0.68 h −1 based on two oxidized humic pools (Roden et al., 2010).Field measurements reported minimum electron transfer of 0.8 mmol charge (eq.) m −2 day −1 generating CO 2 at 0.2 mmol m −2 day −1 (Heitmann et al., 2007).This rate was similar to the small production rate of CH 4 at the investigated bog site.
Based on this limited information, we conceptually modeled the reduction and oxidation of humic substances using first order kinetics (Appendix A, Eqs.A133-136).The initial values of the EA (electron acceptors) and ED (electron donors) pools in the humic substances are calculated from the SOM C pool by a ratio of 1.2 eq.(mol C) −1 (Roden et al., 2010).The initial electron accepting capacity used in the model was ca.2000-4000 mmol charge m −2 for the upper 60 cm of peat per m 2 , which is close to the capacity of 2725 mmol charge m −2 derived from a drying and rewetting experiments in a minerotrophic fen (Knorr and Blodau, 2009).
In the model electron acceptors are renewed via two mechanisms: direct oxidation by O 2 due to WT fluctuation in the only temporarily saturated layers and microbially mediated electric currents through the peat column via an extracellular electron transfer (I nanowire ).While the first mechanism is well documented (Knorr and Blodau, 2009), the second is speculative.It relates to the observation that even in deeper peats, that are not affected by influx of oxygen or other inorganic electron acceptors, CO 2 seems to be net released in excess of methane (Beer and Blodau, 2007).This finding has remained enigmatic because excess CO 2 release would be impossible from a stoichiometric point of view when organic matter with oxidation state close to zero is respired and other, more reduced decomposition products, in particular molecular hydrogen, are not concurrently released.A relevant accumulation of molecular hydrogen has, to our knowledge, not been observed in affected peats.Anaerobic methane oxidation may appear as a way out of the dilemma; however, also this process would depend on the elusive electron acceptor (Smemo and Yavitt, 2011).
Recently an extracellular electron transfer was described that has the potential to solve this enigma.Microorganisms in soils and sediments were first detected extracellularly utilizing electrons from redox active species, such as HS, and Fe (III) (Lovley and Coates, 2000).The term "microbial nanowire" was proposed later for this extracellular electron transfer (Reguera et al., 2005).Recently the process was demonstrated to occur in marine sediments over macroscopic distances (Nielsen et al., 2010).The authors suggested that electrons can extracellularly flow in interconnected networks of "nanowires" so that oxidation and reduction process are spatially separated from each other.In our case the oxidation process releasing CO 2 would proceed deeper into the peat, whereas the reduction reaction would take place near the peatland surface where oxygen is present.We suppose that this mechanism may be the reason for some of the frequently observed CO 2 production that is unrelated to physical supply of an electron acceptor deeper into the peat.Not knowing about mechanistic detail in peats, we conceptualized this process by simply calculating an extracellular electron current in the peat and using Ohm's law for the anoxic layers (Appendix A, Eq.A137).Peat electron flow resistance (R) is determined by inverse modeling based on the resistance constant definition and corrected for soil moisture under the assumption that air-filled pore space cannot conduct electrons (Appendix A, Eq.A142).The parameter n peat ( Ù m) is the specific resistance of the material and l is the layer depth (m).Electron current in mA was then converted to mmol by the Avogadro constant (NA) and the Faraday constant (F) (96 490 Coulombs mol −1 ) (Appendix A, Eq.A137).To make this process work, electrochemical potential gradients (dEh) that drive the flow between adjacent layers are needed.In absence of meaningful measurements of redox potential of peat we calculated such a gradient from a measured redox potential gradient in the Mer Bleue Bog that was given by concentration depth profiles of dissolved H 2 , CO 2 , and CH 4 .We assumed that the redox potential gradient of this redox couple represents the minimum depth gradient in electrochemical potentials being present.Using the Nernst equation for the reaction 4H 2 (aq) + CO 2 (aq) → 2H 2 O (l) + CH 4 (aq) (Appendix A, Eqs.A138-A141), concentration profiles were converted into electrochemical potential gradients with depth.H 2 concentration was measured by Beer and Blodau for the Mer Bleue bog (2007) (Table S4).
In the model the electron flow through the peat towards the peatland surface is used to reoxidize H 2 S to sulfate and DOM-QH 2 to DOM-Q at larger depths.These species are the reduced again, producing the needed "excess" CO 2 in the process and lowering rates of methanogenesis, respectively (Appendix A, Eq.A136).The rate constant of sulfate reduction was adjusted to the suggested range of the SO 2− 4 reduction rates based on the S deposition on the site at 0.89 mmol S m −3 day −1 (Vile, 2003a).The same thermodynamic inhibition concept as used to model methanogenesis was applied also to bacterial sulfate reduction (Appendix A, Eq.A129).
Both CO 2 and CH 4 are in equilibrium between gaseous phase and dissolved phase obeying Henry's Law (Appendix A, Eqs.A100-A103).The efflux of C and N are through runoff and advection in dissolved phase and in gaseous phase from the soil surface.Diffusion follows Fick's law with moisture-corrected coefficients in the saturated layers and was modeled as step functions in the unsaturated layers where diffusion accelerates by orders of magnitude for gases (Appendix A, Eqs.A104-A107).CH 4 also escapes from the soil via ebullition and plant-mediated transportation (Appendix A, Eqs.A115-A120).Ebullition occurs in saturated layers once CH 4 level exceeds the maximum concentration CH 4,max .The parameter CH 4,max is sensitive to temperature and pressure (Davie et al., 2004), with a base maximum CH 4 concentration at 500 uM, which is the value for a vegetated site at 10 • C in Walter et al. (2001).The ebullition of CH 4 releases the gas to the atmosphere without it passing through the unsaturated zone.In the rooted layers, graminoids transport CH 4 at rates that are determined by the biomass of the graminoid roots.A percentage of 50 % of the CH 4 are oxidized to CO 2 during the plantmediated transportation by the O 2 in plant tissues (Walter et al., 2001).The CH 4 oxidation in the oxic layers was modeled using temperature-sensitive double Michaelis-Menten functions (Segers and Leffelaar, 2001) ) in the oxic and anoxic layers, respectively.During nitrification, the fraction of N loss as NO (rNO nitri ) is 0.1 %-4 % day −1 with a mean value of 2 % (Baumgärtner and Conrad, 1992;Parsons et al., 1996).For N 2 O (rN 2 O nitri ) this value is smaller at 0.1 %-0.2 % day −1 (Ingwersen et al., 1999;Breuer et al., 2002;Khalil et al., 2004a).We used similar values as in the model DNDC for acid ecosystems, where rN 2 O nitri was 0.06 % and rNO nitri was 0.25 % (Li and Aber, 2000).Both nitrification and denitrification are regulated by temperature, moisture, and pH.Moisture is the dominant control for nitrification and an effective control for denitrification (Linn and Doran, 1984;Riedo et al., 1998).
In an acidic environment, nitrification was detected to cease below pH of 4 and reached a maximum at a pH of 6 (Lång et al., 1993).The optimal range of pH for denitrification was suggested to be from 6 to 8 (Heinen, 2006).Temperature factors were empirically modeled, using the equation in DNDC (Li and Aber, 2000) for nitrification and the common formalism equation in NEMIS (Johnsson et al., 1987;Hénault and Germon, 2008) for denitrification.
3 Model application

Site description
The model was applied on the Mer Bleue (MB) Bog for a period of 6 yr from 1999 to 2004 to evaluate the simulation performances of WT dynamics, carbon fluxes, soil water DIC and CH 4 concentrations and C and N budgets against observations.
The Mer Bleue Bog (45.51 • N, 75.48 • W) is a raised acidic ombrotrophic bog of 28 km 2 located 10 km east of Ottawa, Ontario.The bog was formed 8400 yr ago as a fen and developed into a bog between 7100 and 6800 yr BP.The peat depth varies from 5 to 6 m at the center to < 0.3 m at the margin (Roulet et al., 2007).The vegetation coverage is dominated by mosses (e.g.Sphagnum capillifolium, S. angustifolium, S. magellanicum and Polytrichum strictum) and evergreen shrubs (e.g.Ledum groenlandicum, Chamaedaphne calyculata).Some deciduous shrubs (Vaccinium myrtilloides), sedges (Eriphorum Vaginatum), black spruce (Picea marinana) and larch also appear in some areas (Moore et al., 2002).The annual mean air temperature record from the local meteorology station is 5.8 degrees and the mean precipitation is 910 mm (1961-1990 average;Environmental Canada).The coldest month is January (−10.8 • C) and the warmest month is July (20.8 • C) (Lafleur, 2003).

Application data and initialization
Inputs required are geographic location and local slope of the site, daily precipitation and PAR, daily snow-depth record, annual average and range of air temperature, atmospheric CO 2 , CH 4 and O 2 levels, annual N load and vegetation type of the site (Table 2).
Observed C fluxes, water table depth, and the depth profiles of temperature and moisture with 5 s to 30 min intervals were obtained from fluxnet Canada (http://fluxnet.ccrp.ec.gc.ca) and averaged to daily values.Fluxes were determined using micrometeorological techniques, and gaps shorter than 2 h were filled by linear interpolation between the nearest measured data points.Longer gaps were filled by repeating the corresponding period of time from the closest available dates.Other data sets for model evaluation were obtained from a range of the published literature.The spin-up (initiation) of the model was conducted with initial values obtained from literature (Table S4) and the meteorological and geophysical boundary conditions (Table 2) from 1999 to 2004 obtained from fluxnet Canada.The time series was repeated every 6 yr until the model approached its steady state after a period of longer than 100 yr.The obtained values of state variables were used for the actual model application and evaluation.Most parameters were obtained from literature for bogs or peatlands in general, or calibrated for the ranges from measurements, or in line with the values used in previously published models.In total, 29 out of 140 parameters were calibrated and ranked from 3 to 1 based on their origin and descending confidence in their accuracy and correctness (Tables 3 and 4).Parameters in category 3 were calibrated with comparison to similar parameters in references; parameters in category 2 were calibrated in comparison to conceptually related parameters in references; parameters in category 1 were unavailable in literature and thus were calibrated without references (Table 4).

Results
We ran the parameterized, initiated model for 6 yr from 1999 to 2004 and evaluated the simulation results of WT depth, and depth profiles of soil temperature, moisture and O 2 to assess the ability of the model to generate environmental controls on C and N cycling.The simulated C and N pool sizes, transfer rates and fluxes were compared with six years of continuous measurements to evaluate the capability of the model in quantifying C and N pools and cycling rates.We also conducted sensitivity analysis for the key factors (e.g.temperature, precipitation, N deposition) and a range of uncertain calibrated parameters (e.g.potential decomposition rate of the soil organic matter).This demonstrated the sensitivity of the model to N availability and climate controls, which shows the potential for applying the model to longterm N fertilization and N deposition and climate change studies.As statistics for evaluation we chose the root mean square error (RMSE), linear regression coefficient (r 2 ), and the index of agreement (d) (Willmott, 1982).

WT depth, soil temperature and moisture
Simulated daily average soil temperature was plotted against measured temperatures in hummocks at 0.05 m and 0.8 m depth (Fig. 5a).The simulations agreed well with the observations and showed degrees of agreement (d) of 0.97 and 0.95, and RMSE of 3.23 and 1.70 degrees, respectively.However, the model failed to simulate the observed deviation from the sinusoidal temperature curve when snow was not present in the winter of 2003, implying other controls on soil temperature that are currently missing in the model.In general, the simulated WT depth showed good agreement with the observed data, with a degree of agreement (d) of 0.98 and RMSE of 0.06 m (Fig. 5b).The largest deviation was from mid-July to early August of 1999, when the simulated WT depth for some days reached the maximum depth and was more than 20 cm below the observed WT depth.From 1999 to 2002, WT depth elevation was underestimated during seasonal changes from summer to fall when the deviations of more than 10 cm occurred for 10 to 30 days.These disparities were likely owed to the simple bucket model structure that lacks processes of water transfer that buffer variations in water content.
Considering the large variation of soil moisture between hummocks and hollows, we compared the simulation at 0.2 m and 0.4 m depth with the observations in hummock and hollows, respectively (Fig. 5c).The seasonal dynamics were well captured and the 0.4 m simulation agrees with the observation strongly.However, the simulated volumetric water content at 0.2 m was systematically overestimated by 0.1 to 0.2 in summers and up to 0.5 for the wettest year in winter.
Large spatial in situ variability of observed volumetric water content might be one of the reasons for this large discrepancy, as the simulated values are similar to other measurements in hummocks in the Mer Bleue Bog during even drier years (Wendel et al., 2011).

Daily Carbon fluxes
Gross ecosystem production (GEP) was calculated as the sum of simulated gross primary production (GPP) of all PFTs (Fig. 6a).The simulated ecosystem respiration (ER) was the release of CO 2 gas from the peat surface, which included autotrophic respiration (AR) in shoots and roots of plants and the heterotrophic respiration (HR) of microorganisms in the soil (Fig. 6b).Net ecosystem exchange (NEE) was calculated as the difference between ER and GPP (Fig. 6c).
Overall, the simulated GPP, ER and NEE captured the seasonal dynamics and the magnitudes of the C fluxes.The maximum simulated daily GPP was 5.96 gC m −2 day −1 and occurred in the driest year 1999, which is similar to the maximum observed 6.80 gC m −2 day −1 .The simulated starting dates of spring PSN ranged from day 79 (2000) to day 99 (2001), with an average date of day 90.These values fell in the reported range from day 86 to day 101 (Moore et al., 2006).The simulated starting dates of PSN in 2001 and 2003 were at day 99 and 84, which was two days earlier than in field observations.The average difference between simulated and observed GPP was 0.43 gC m −2 day −1 , which was slightly larger than the calculated mean error of GPP (±0.11 gCO 2 m −2 day −1 ) in measurements (Moore et al., 2006).Statistical analysis revealed a root mean square error (RMSE) of 0.73 gC m −2 day −1 and a degree of agreement (d) of 0.95 (Fig. 7a).However, there were a few days when the simulation errors were large, among which the maximum underestimation was 3.68 gC m −2 day −1 on 31 July in 2000 and the maximum overestimation was 3.21 gC m −2 day −1 on 23 May 2002.
Daily CH 4 flux was simulated from 1999 to 2009 in order to compare with the observations from 2004 to 2008.Simulated daily CH 4 flux covered a wide range from 0 to ca. 170 mg m −2 day −1 .Seasonal patterns were stronger in wet years, such as 2004 and 2006, when the fluxes reached a maximum in mid-summer.In the dry years (e.g. 2005, 2008), summer peaks were lacking and the maximum fluxes occurred during one day in late spring and early summer due to degassing when the water table quickly declined (Fig. 8a  and b).The instantaneous degassing in the model was caused by the release of CH 4 stored in each 5 cm layer that entered the unsaturated zone.Subsequently the CH 4 fluxes fell to very small values due to limited production and increased methane oxidation during summer.The simulated CH 4 flux agreed with the observed range from April to mid-May and was underestimated in summer (Fig. 8b).

Dissolved CH 4 CO 2 and O 2 concentration
The simulated daily concentration of dissolved CH 4 and CO 2 was plotted against depth for 2002 to evaluate the model output of belowground respiration (Fig. 9a  and b).Both dissolved CH 4 and CO 2 accumulated with depth and showed clear seasonal dynamics with the seasonal WT fluctuation.Concentration of dissolved CH 4 increased from <0.1 mmol L −1 around the WT at 0.35 cm to ca. 0.6 mmol L −1 at 80 cm depth in January and to ca. 0.5 mmol L −1 at 90 cm in October.Concentration of dissolved CO 2 increased from <0.1 mmol L −1 around the WT to ca. 3.5 mmol L −1 at 70 cm depth in January and to over 6 mmol L −1 in October.The maximum concentration in deep layers was ca.7 mmol L −1 dissolved CO 2 and 0.6 mmol L −1 dissolved CH 4 , respectively, close to the observed ranges (Beer and Blodau, 2007;Beer et al., 2008).Figure 9c illustrates the profile of dissolved O 2 concentration for the year 2002.The dissolved O 2 was depleted rapidly below the WT, where concentration decreased from ca. 0.3 mmol L −1 at around the WT in January to ca. 0.1 mmol L −1 in October.Summer O 2 concentration around the WT was lower than the rest of the year, due to the alteration of Henry's law constant of O 2 by the increased summer temperature.Oxygen in soil was consumed by two processes in the model: organic C oxidation and methane oxidation.The annual consumption of O 2 in methane oxidation was between 5 % and 7 % of the annual input of O 2 from the atmosphere that diffused into the soil during the simulation period.Therefore methane oxidation was not an insignificant sink of oxygen, yet it was not highly important either.

Annual C budget
We calculated an annual C budget (Fig. 10a) based on the 6 yr mean of annual simulated pool and flow rates (Table S1).Annual GPP ranged from 513 gC m −2 yr −1 in the second wettest year 2000 to 609 gC m −2 yr −1 in one of the dry years 2001.Similar to the 550 gC m −2 yr −1 of GPP in the conceptual C budget model for the Mer Bleue Bog (Moore et al., 2002) Roulet et al., 2007), which was based on 8 yr of observations from 1999 onwards.The model simulated an annual CH 4 emission of 4 gC m −2 yr −1 of which 83 % stemmed from graminoid mediated emission.Emission of CH 4 during the wet years of 2002 and 2004 were higher than in the dry years, as is the general trend observed in the Mer Bleue Bog and in other peatlands (Roulet et al., 2007).The simulated DOC export was 15 gC m −2 yr −1 , which was in agreement with the estimated 14.9 (±3.1) gC m −2 yr −1 from 5 yr of runoff and 3 yr of DOC concentration measurements at the site.

Annual N budget
An annual N budget for the Mer Bleue Bog is illustrated based on the 6 yr average of simulated values (Fig. 10b, Table S2).The wet annual N deposited from the atmosphere was 0.81 gN m −2 yr −1 onto the peatland.About 95 % of the deposited N was absorbed by mosses right away.Nitrogen in the plants was associated with the plant biomass and composition, which both changed little over the 6 yr.Annually, mosses exported 0.82 gN m −2 yr −1 in litter and 0.02 gN m −2 yr −1 in exudates to the soil N pools.For vascular plants these fluxes were 2.97 gN m −2 yr −1 and 0.02 gN m −2 yr −1 , respectively.N uptake was 1.68 gN m −2 yr −1 , mostly by shrubs as NH + 4 , and only 0.3 % of N uptake occurred in form of DON.N 2 fixation was 0.96 gN m −2 yr −1 .Considering N uptake, N litterfall and N exudation, vegetation thus lost 0.38 gN m −2 yr −1 , which represents 2.5 % per year over the simulation period.The NH + 4 pool was smaller than the annual production and uptake, implying a fast turnover of NH + 4 in the soil.Other dissolved N pools (NO − 3 , N 2 O and NO) were 3 to 8 magnitudes smaller than the NH + 4 pool in the model, and N 2 O emission was negligible.Export of DON and DIN through water runoff was also very small and occurred at rates of 0.04 gN m −2 yr −1 and 0.01 gN m −2 yr −1 , respectively.Overall, the OM pools received ca.3.83 gN from plant litter production and exudation and lost 1.91 gN and 0.05 gN by mineralization and runoff annually, which lead to an overall accumulation of 1.43 gN m −2 yr −1 in the peat.

Sensitivity analysis
Sensitivity analysis is useful in quantifying the model responses to changes in environmental drivers and other parameters.We ran a series of simulations by adjusting key environmental variables, such as precipitation, air temperature and N deposition.Variations of these parameters were chosen to be within the possible range of variability in temperateboreal peatland ecosystems.We also adjusted parameters that are most uncertain and potentially influence C and N cycling in peatlands, such as Q 10 values and the rate constants of the decomposition of SOM.The sensitivity of key C and N fluxes, pools, and cycling rates, including GPP, AR, ER, HR, NEE, NECB, and C and N sequestration rates in the soil organic matter, were examined.The sensitivity was tested by imposing changes in air temperature between −1 and +5 with increments of 2 • C, and changes in precipitation between −30 % and +30 % with increments of 15 %, which were in line with the scenario predictions of future climate (IPCC, 2007).The sensitivity to N deposition that covered the N deposition range in Europe was tested by imposing N input at 0.2, 1.4, 2, 2.5 and 3.2 gN m −2 yr −1 .The sensitivity to Q 10 of labile (Q 10,L ) and recalcitrant (Q 10,R ) soil organic matter were tested for −40 % and +40 % of the ambient value, respectively.The potential decomposition constant k was tested with −25 % and +25 % of the ambient k for labile (kpot L ) and recalcitrant (kpot L ) in the sensitivity tests.

Y. Wu and C. Blodau: PEATBOG
The simulations were conducted for six years and averaged to compare with the baseline simulations (Table 5).
The sensitivity analysis showed that heterotrophic respiration was the most sensitive process in C cycling with regard to air temperature.Temperature increase had a negative effect on the production of moss and a positive effect on the production of vascular plants, suggesting a favoring of vascular species in a warmer environment.The increase in AR in vascular plants with increasing T was greater than the increase in production of vascular plants, which led to a negative effect on NPP.In the model, the Q 10 of respiration in plants was smaller than the Q 10 of photosynthesis, suggesting that other controls constrain primary production apart from temperature, such as N availability and soil moisture.The sensitivity of HR to temperature was greater than that of AR, resulting in preferential C loss from peat rather than from plant respiration with increasing temperature.The impact of temperature on ER was larger than on GPP and entailed a higher sensitivity of NEE to temperature as well.Although less CH 4 , DOC and DIC was exported when temperature was increased, NECB declined due to the greater change in NEE.Carbon sequestration was very sensitive to temperature in the model, and an increase of 1 degree in air temperature would turn the modeled peatland from a C sink into a C source.Nitrogen sequestration was also negatively affected by temperature, but to a lesser extent.
The processes GPP, AR, HR were less sensitive to precipitation than to temperature.This was not the case for the export of dissolved C and CH 4 fluxes.Decreasing precipitation promoted primary production and autotrophic respiration in vascular plants, while inhibiting the production of mosses.Increasing precipitation more strongly raised NPP in shrubs than in mosses and had a negative effect on graminoids, suggesting vice versa that graminoids were more tolerant to dryness than shrubs and mosses.The increased NPP in shrubs resulted mostly from changing respiration rather than from gross primary production.Respiration in the model has a stronger dependency on soil moisture than GPP.In the analyses, HR was more sensitive to temperature and precipitation than AR and NPP, and it was more sensitive to temperature than to precipitation (Table 5).A decrease in precipitation by 30 %, corresponding to a decline of annual mean WT depth by 7 cm, led to an HR increase of 11 %.In contrast DIC and DOC export declined by 36 % and 66 %, respectively.The decrease of dissolved C exports was owed to the diminished runoff at lower WT position, despite more production of dissolved C with raised HR.As expected, CH 4 flux was strongly positively related to precipitation.In contrast, elevated temperature decreased CH 4 emission in the model through the lowered WT depth (Table 5).
Interestingly, the sequestration rate of C was similarly sensitive to precipitation and to temperature, while the N sequestration rate was much more sensitive to precipitation than to temperature.A decrease in precipitation by 30 % caused a decrease in C sequestration rate by 19 %, which is compara-ble to the effect of an increase in temperature by 3 degrees.Meanwhile, the N sequestration rate decreased by 46 % with the change in precipitation and by 10 % with the change in temperature.This outcome resulted from the different mechanisms by which precipitation and temperature control the decomposition of soil organic matter.In the model, lowering the WT position via precipitation stimulated the decomposition rate of labile and recalcitrant soil equally.On the other hand, the temperature increases primarily the decomposition of recalcitrant OM due to a larger decomposition Q 10 of this pool.As recalcitrant soil is present mostly in the deeper layers and contains less N, the temperature effect on N sequestration was weakened.Therefore, if recalcitrant SOM is more sensitive to temperature than labile SOM, as suggested by many (Davidson and Janssens, 2006;Conant et al., 2008;Craine et al., 2010;Karhu et al., 2010), the function of peatlands as N sinks will be more impaired than in predictions on models with equal Q 10 values for labile and recalcitrant SOM.
Nitrogen deposition levels affect mostly plant related C fluxes rather than soil derived fluxes.The sensitivity of GPP to N deposition was greater than to precipitation and temperature.Overall, the model suggests a strong promotion of graminoids over shrubs and mosses when the N deposition increases.The effect of N on both GPP and NPP was stronger in graminoids than in shrubs and mosses, due to the different N use strategy of the PFTs in the model (Table 5).Graminoids have advantages because faster turnover rates allow for instantaneous response to changes in N availability in the plant-soil system.In comparison, shrubs and mosses cycle N in a more conservative manner and need lower levels of N to keep photosynthesizing, hence these PFTs react more slowly to increases in N availability.The NPP of graminoids increased non-linearly with the N deposition level, by 70 % with a 150 % increase and 560 % by a 300 % increase in annual N deposition (Table S3).This finding implies other constraints on the NPP of graminoids at low N deposition levels.The main constraint was very likely N filtration by mosses, which was alleviated when mosses became N saturated at higher N deposition levels.
The NPP of shrubs was highest at moderate N deposition level of 2.6 gN m −2 yr −1 , probably due to increased shading effects from the faster expansion of graminoids with more N deposition (Table S3).The NPP of mosses was negatively affected by N deposition, and only a slight promotion of GPP occurred when N deposition was slightly raised.Very different from the effects of the climatic drivers, N deposition levels had hardly an effect on HR.Other C effluxes, including dissolved C export, CH 4 flux and AR were also less sensitive to N deposition than to temperature and precipitation.As GPP and ER were both positively affected by increasing N, the NEE, NECB and C sequestration rate of peat were not very sensitive to N deposition.In contrast, N sequestration in soil organic matter showed a strong positive relation to N deposition level.(Jørgensen and Bendoricchio, 2001).Annual C fluxes (unit: gC m −2 yr −1 ) averaged over 6 yr from 1999 to 2004 were compared per change of air temperature (unit: • C), precipitation (unit: m day −1 ), N deposition level (gN m −2 yr −1 ), Q 10 (no unit) and kpot (potential decomposition constant, unit: day −1 ) of the labile and recalcitrant Peat.

Parameters
Air Temperature Precipitation N deposition Processes in the model were generally more sensitive to changes in parameters related to the recalcitrant OM fractions (Table 5).Plant-derived C fluxes were little sensitive to Q 10,L , Q 10,R and kpot L , but moderately sensitive to kpot L .The effects of kpot L on GPP occur through changes in N availability in the peat, which varies according to the decomposition rate of the recalcitrant soil.The processes HR, NEE, NECB and the sequestration rates of C and N in soil showed greater and significant sensitivity to kpot L and Q 10,R than to kpot L and Q 10,L , showing the importance of the recalcitrant SOM pool for HR.In the short term, the process most sensitive to all varied factors other than kpot L was the net ecosystem carbon balance (NECB).

Nitrogen saturation
Increased N deposition has been observed to change vegetation composition and the C and N retention in mosses, vascular species, and peat (Lamers et al., 2001;Xing et al., 2010;Bragazza et al., 2012).The model was in part designed for quantifying changes in PFTs and for identifying the threshold of N deposition level where N saturation occurs in mosses.
To study the plausibility of the model behavior we carried out a 40 yr simulation with raised atmospheric N input (Fig. 11).We adjusted the N deposition to 1.5 gN m −2 yr −1 , which is the intermediate N deposition in the sensitivity analysis and has been suggested to be the critical load of N for mosses (Vitt et al., 2003).The C and N pools in PFTs showed a delay in responses to elevated N deposition (Fig. 11a and  b).The fraction of deposited N absorbed by mosses remained steady for the first 12 yr until the N content reached 0.02 gN g −1 biomass (Fig. 11d).Above this content level, the fraction of N retained by mosses declined rapidly and excess N entered the pore water.As a result, only then did the fraction of deposited N retained in vascular plants and peat increase and peaked after ca.20 yr (Fig. 11c).
Nitrogen mineralization rates increased immediately after raising N deposition, because of the elevated litter production in plants and exudation of mosses (Fig. 11f).Output of N from the model ecosystem was about 5 % of the total N input from deposition and N 2 fixation, and was continuously increasing after moss filtration of N became less effective (Fig. 11f).
One of the important findings of this exercise was that total biomass and total NPP remained comparatively stable, while the plant composition of biomass and NPP changed greatly (Fig. 11a and e).The moss cover was completely diminished, while graminoids started to expand with higher N availability in the soil water and eventually became the dominant PFT.An increase in the labile fraction of SOM was a further consequence because invading vascular plants produce more labile litter in the model.Owing to both the increased litter inputs from the vegetation and raised litter decomposability, the sequestration rate of C in soil first accelerated but then slowed after the NPP had peaked (Fig. 11e).

Carbon fluxes and environmental controls
The fluxes GPP, ER and NEE are the essential components in C cycling that express the ability of peatland ecosystems in assimilating and dissimilating C and exchange the element with the atmosphere.Overall, the model simulations showed good agreement in daily C fluxes, belowground C concentration and annual C and N budgets with empirical data.However, a bias occurred towards underestimating simulated GPP (i.e.slope = 0.936), underestimating simulated ER (i.e.slope = 0.806) and overestimating simulated NEE (i.e.slope = 1.166).These biases are within the bias range of the other models that primarily focus on C cycling (e.g.MWM, PCARS).The model performance differed in that in MWM and PCARS the simulated ER was overestimated, while it was underestimated in the PEATBOG model.
The 6 yr averaged annual GPP demonstrates the ability of the model in simulating overall productivity, as only a small deviation of 5 gC m −2 was recorded against an empirically determined large average GPP of 550 gC m −2 at the site (Moore et al., 2002).Also the trends in interannual variation of GPP with precipitation and temperature were largely met.Noteworthy is for example the decline in GPP in the extremely dry year 1999, when dryness had a large impact on the GPP of mosses, and the high GPP in the warm and wet year of 2001 (Figs.5b and 6a).While overall model performance was good some deviation from empirical measurements were illustrated by the analysis as well.Annual GPP was overestimated by 32 to 85 gC m −2 yr −1 from year 2000 to 2003 and underestimated by 70 to 123 gC m −2 yr −1 for the remaining years by the model simulations (Table 6).The discrepancy of annual GPP simulations ranged from 7 % to 18 % and was not significant (P = 0.737, n = 2192).The simulated GPP fraction of shrubs was 70 %, ranging from 66 % in the simulated wettest year of 2004 and 78 % in the driest year 1999.This range was similar to the model output of MWM that ranged from 61 % to 67 % (St-Hilaire et al., 2010) and smaller than the shrub related fraction of GPP of 80 % to 85 % reported from the PCARS model (Frolking et al., 2002).Inter-annual variation of GPP for PFTs was corroborated by observation (Bubier et al., 2003): GPP of mosses increased from dry to wet years from 4 % to 48 %, whereas GPP of shrubs was at its lowest levels in the wet years.In comparison to other models (St-Hilaire et al., 2010;Dimitrov et al., 2011), the inhibition of GPP of shrubs due to dryness is less effective in our model.
On the daily timescale some weakness of the model in responding to weather conditions became visible.In general, the simulated GPP was deficient in capturing short-term extreme fluxes.All large underestimates (>2 gC m −2 day −1 ) in the GPP simulation occurred during mid-summer in the two wet years 2000 and 2005, when GPP in the peatland was larger than 5 gC m −2 day −1 , except for two days in late summer.The likely reason for the lack of adequate model performance during this time are the maximum photosynthesis rates that are set for each PFT and the impossibility to cover the daily observed extreme values that were averaged from half hour records in the measurements.This disadvantage also occurred in other models with maximum rate settings that are based on the Farquhar photosynthesis model (e.g.MWM).We also noticed that most of the underestimates that occurred in 2004 were associated with frequent heavy precipitation that raised production instantly.In the model, the production of mosses is the only PFT that reacts to precipitation directly through the water content in the capitulum of mosses.The indirect controls of precipitation on the production of vascular plants via WT depth is likely the reason of the underestimated promotion of photosynthesis by frequent precipitation, especially when other controlling factors (i.e.temperature, light) are within the optimal range.For example, a peak of measured daily GPP was observed during late July 2004, during one of the periods that underestimated GPP.At this time precipitation was continuous at >10 mm day −1 and temperature was within an optimal range (20 ± 3 • C).
The overestimation of GPP mainly occurred during late May to early June in the dry years (2001 to 2003) when PAR was comparably strong (>600 umol m −2 s −1 ).During those days, the model predicted GPP of mosses and shrubs to reach a level above 1.2 gC m −2 day −1 and 2 gC m −2 day −1 , respectively.Daily measured GPP in the Mer Bleue Bog was found to be significantly albeit weakly related to PAR (P < 0.001; r 2 = 0.19) (Moore et al., 2006).In the this relationship is significantly stronger (r 2 = 0.75), due to neglecting the non-linearity of leaf response to light in the integration of canopy photosynthesis using just Beer's law.The non-linearity of leaf response to light is related to the diurnal effects on the canopy.It includes for example optimized nitrogen distribution in plant canopies, different responses to light in sun and shade leaves, and variation of stomatal conductance with light levels (Thornley, 2002;Hikosaka, 2003).Late May to early June was also the period when new biomass is built up, which affects the distribution of N within the plants.For example, both total N content and chlorophyll-a concentration in evergreen shrub foliage were low in spring and increased steadily to early June, as shown in measurements (Moore et al., 2006).The model lacks separated N pools in foliage and stems, where N content could show great variations due to phenology, which might be the reason of the overestimation of GPP in late spring.
The fluxes ER and NEE represent the gross and net release of CO 2 from peatlands, and largely determine the C balance of the ecosystem.The model reproduced the composition of ER, where HR contributed half of the total ER, while the other half was almost equally shared by AR in shoots and AR in roots, as approximately suggested by field measurements (Moore et al., 2002).However, the standard deviation of the simulated annual ER and NEE was larger than that in field estimates (50 % and 40 %), suggesting a larger inter-annual variation than measured in the field.The modeled annual ER ranged from 430 gC m −2 yr −1 to 573 gC m −2 yr −1 , with an average of 526 gC m −2 yr −1 (Table S1), which is close to the flux quantified as 461 gC m −2 yr −1 (Lafleur et al., 2001).The annual discrepancy ranged from 3 % to 17 %, with an exception of 25 % in 2004, when the highest summer WT occurred (Table 6).The underestimation of ER was probably caused by the simulated WT depth (Fig. 5b) that was 5 to 10 cm higher than measured in summer when both autotrophic respiration (AR) and heterotrophic respiration (HR) were potentially high.The modeled NEE showed similar interannual patterns to ER with the annual error being between 35 gC m −2 yr −1 to 18 gC m −2 yr −1 .The largest deviation of simulated NEE from measurements was 106 gC m −2 yr −1 in 1999, when GPP was under-and ER overestimated.
The ER was also overestimated from 1999 and 2003 (Fig. 6b).To identify the reasons, we calculated the deviation between measured and model daily ER, and regressed it against the deviation of measured and modeled daily temperature.According to this approach the overestimate of ER from 1999 to 2003 could be explained by an overestimate of soil temperature (r 2 = 0.26), especially during summer (r 2 = 0.68).Both ER and HR were strongly correlated to soil temperature at 0.2 cm depth with r 2 of 0.88 and 0.83, respectively (n = 2193).The strong temperature dependence of ER and HR was associated with the Q 10 values used in the model for the temperature effects on HR rates.Different from other models, where Q 10 values were set to 2 for microbial respiration in soil, Q 10 value for the decomposition of recalcitrant OM (3.3) was set to be larger than for labile OM (2.3).These Q 10 values were in line with some of the most recent results (Davidson and Janssens, 2006;Conant et al., 2008Conant et al., , 2010;;Karhu et al., 2010), their application implies a stronger increase in C loss from peatlands in a warmer climate.It has to be noted that some have assumed the value of Q 10 for labile OM to be larger (Liski et al., 1999;Giardina and Ryan, 2000;Thornley and Cannell, 2001) than or similar to (Fang et al., 2005) that of recalcitrant OM; in this case climate change effect on NEE may not be as extraordinary as has been anticipated otherwise.The sensitivity analysis on Q 10 and potential decomposition rates for our model highlighted the importance of the recalcitrant OM over the labile OM for the C cycling in peatlands (Table 5).
The Q 10 values derived from the first order exponential equations of the simulated ER and HR were only 2.56 and 1.97, respectively.The Q 10 for HR was thus smaller than either of the Q 10 for labile or recalcitrant OM, revealing the importance of other factors that confound the temperature response of HR.The WT depth was the most important factor affecting the calculated Q 10 values with r 2 of 0.75 between the average summer WT depth and the calculated Q 10 values.In summer, the low soil moisture in the most upper peat layers counteracted the potential enhancement of respiration by temperature.Nevertheless, the sensitivity analysis suggested a lesser effect of WT depth than of soil temperature on CO 2 fluxes (Table 5).The daily simulated WT depth moderately correlated with ER (r 2 = 0.51, n = 2192), with r 2 ranging from 0.19 in the wet year 2000 to 0.79 in the dry years.Although stronger than reported from empirical studies, this relationship was in a broad agreement with field results as far as the trend of tighter correlation in dry years goes (Bubier, 2003;Lafleur et al., 2005b;Blodau et al., 2007).
The CH 4 fluxes modeled with our novel thermodynamickinetic approach were in a reasonable range but smaller and their seasonal pattern less pronounced than obtained with chamber measurements at the Mer Bleue Bog (Moore et al., 2011).We attribute this difference to the variability of in situ plant cover and a higher mean water table position of the 12 gas flux collars of the field study.The collars were not only located in hummocks and lawns but also in hollows.The observed average WT depth was −35 (±8.4 cm) for the 12 collars from 2004 to 2008, whereas the simulated average WT depth was −41 cm for the same time period.Due to the generally observed exponential increase in emissions with raising water table (Moore et al., 1998), even a small number of sampled wet locations may lead to much larger emissions than simulated in the model, which represents a hummock situation.The large discrepancy after summer was very likely due to the effects of vegetation transport on CH 4 flux, which was the most important control on the CH 4 flux from September to November over 5 yr (Moore et al., 2011).In the model, graminoid cover was less than 1 % during the simulation period, whereas the graminoid cover ranged from 0 to 100 % in the 12 collars.Comparing model output to one of the gas flux collars similar in water table and graminoid cover (collar 8, Table 2, Moore et al., 2011) with daily CH 4 flux between 10 to 100 mg m −2 day −1 , a closer model fit was obtained.In this collar, as in our model, CH 4 emission increased less in summer than in the more grass-rich collars (Moore et al., 2011).
The growing season log 10 values of both daily and annual CH 4 fluxes showed moderately strong relations with WT depth (r 2 = 0.56, n = 2119 and r 2 = 0.45, n = 11) (Fig. 8c  and d).The outliners were the degassing events, which occurred when water table was crossing the boundaries of peat layers in the model.The WT depth during the growing season showed differing effects on CH 4 fluxes in dry years and wet years, as was also found in the field (Moore et al., 2011).According to the model results, the lowest dependence of CH 4 flux on the WT depth occurred in the dry years and the highest dependence in the wetter years.This finding is in conflict with relations obtained from field data, where CH 4 emissions were less related to summer WT depth in wetter years.The annual variation in CH 4 production is less pronounced than in CH 4 fluxes (Table S1); this implies that changes in the transport mode of CH 4 might offset the well-known WT control on methanogenesis.For example, the root biomass of graminoids, that provide conduits for CH 4 transport, was negatively correlated with WT depth and CH 4 fluxes.In the dry years, graminoid root biomass increased with declining WT in the model, due to more reallocation of newly produced biomass to roots for accessing soil water.This adaptation also increased the transport of CH 4 from the deeper peat.Overall, the model was able to simulate the variation of CH 4 fluxes with the change of environmental controls, and revealed some interesting dynamic interactions with ecosystem structure that warrant further analysis in the future.

N budget and N saturation
The simulated N budget identified the Mer Bleue Bog as a currently N-limited ecosystem and sink for the element.The immobilization of deposited N by mosses was at a maximum level of 95 %, including both the retention of N in the capitulum of Sphagnum mosses and indirect retention via their stems.In the simulation of N saturation, the model was able to track the effect of N deposition in different compartments of the ecosystem.The N content in mosses peaked at 0.02 gN g −1 biomass, similar to the field observations of 0.015 to 0.024 gN g −1 biomass (Heijmans et al., 2001;Granath et al., 2009;Xing et al., 2010).The simulated increase in soil organic matter mineralization was in agreement with most fertilization experiments (Bragazza et al., 2006;Breeuwer et al., 2008).It was closely related to a change in peat chemistry, such as reflected in the size of the labile OM fraction in peat and its C/N ratio, as observed in a 7-year fertilization experiment (Bragazza et al., 2012).The model also successfully simulated the maintenance of total PFT biomass and production with dramatic changes in the PFT composition, as observed in many N fertilization experiments (Bubier et al., 2007;Juutinen et al., 2010).
Uptake of DON, which has not been considered in peatland models previously, represented a negligible fraction (ca.0.2 %) of the total N uptake by the roots of vascular plants.However, the turnover rate of DON was extremely high, revealing the strong demand and potential uptake of DON by the roots of vascular plants.The fast turnover rate (Kielland et al., 2007) and the large potential uptake of DON (Kahmen et al., 2009) were previously reported from field experiments on boreal forest and three intermediate N available systems, respectively.The model showed that the primary limitation on the uptake of DON was the DON concentration in the soil water, which was also suggested for boreal forest (Kielland et al., 2006) and for Anthoxanthum odoratum in a fertilized experimental site (Sauheitl et al., 2009).Consequently, DON uptake will be more important when there is more bio-available DON in the soil.Although not shown here, the DON uptake accounted for 16 % of the total N uptake of shrubs after 40 yr of N deposition of 1.5 gN m −2 yr −1 in the N saturation simulation.
The nitrogen saturation simulation further showed that the impact of N deposition developed only after a considerable time lag (Fig. 11).Except for mineralization and N output, the C and N pools and fluxes remained stable in the first 12 simulation years until N became saturated in mosses.Only after that point, the N retention in vascular plants and peat increased dramatically and changed the peatland into grass dominated within 8 yr.A delay of 12 yr in the occurrence of effects of N fertilization reveals the importance of accumulated N deposition rather than annual N deposition.

Conclusions
The PEATBOG model has been developed for the purpose of analyzing coupled C and N cycling on a process-level and a daily to multi-year timescale.Our objective was to conceptually consistently integrate vegetation, soil biogeochemistry and soil water dynamics.The model was further designed to be sensitive to changes in N deposition, temperature and precipitation.PEATBOG thus integrates a vegetation submodel comprising three PFTs with a soil and water biogeochemical model providing high spatial and process resolution.It consistently emphasizes mass balance principles and the dynamic interplay of production, consumption and translocation of materials throughout the ecosystem.PEATBOG is able to generate soil physical conditions and plant composition internally and thus requires only a few site specific parameters on geological location, local climate and initial vegetation composition for simulations.The PEATBOG model was effective in reproducing current C and N cycles in a northern peatland with some weaknesses in displaying correct short-term dynamics of C cycling during extreme meteorological periods.It was adequately sensitive to broader changes in climate and N deposition and reproduced a considerable range of empirical findings related to effects of inter annual meteorological variability and N deposition (e.g. the temperature control on soil respiration, change in PFT composition while total C pool and NPP in plants remained robust).
In this paper we presented the components and structure of the model and evaluated the general model performance and sensitivity.The sensitivity analyses and the simulation of increased N deposition demonstrated the model's utility in analyzing the effects of climate change and N deposition on the C and N cycles of peatlands.The analyses further illustrated its usefulness in hypothesis-building that could assist in designing empirical studies examining ecosystem changes over the long term.
In terms of application, the model is suitable for investigating the mechanisms of observed changes in peatland C and N fluxes due to changes in meteorological drivers and N input.Alternatively, the model could be a tool for assessing long-term scenarios of global change.The multi-layer structure of the soil submodel also allows for the integration of other belowground processes in the future, such as SO 2− 4 reduction, to explicitly model CH 4 production on account of the competition among electron acceptors.Although the CH 4 production was modeled conceptually from an electron competition perspective, which we did not detail in this paper, it also produced reasonable annual fluxes and depth profiles of CH 4 concentration, which hold promise for future analyses of CH 4 dynamics.

Figure
Figure 4. Schematic soil water flow

5Figure 5 .Fig. 5 .
Figure 5.Time series of observed daily average (symbols) and daily simula temperature (a),water table depth (b),and volumetric water content (c) for 1 bars in (b) indicate observed daily precipitation records.

7Figure 7 .Fig. 7 .
Figure 7.The scatter plot of observed and simulated daily GPP and ER (a) a 2004, with the best fit relationship (dashed line) and the 1:1 line (solid black li

FigureFig. 8 .
Figure 8.Time series of (a) simulated daily CH 4 flux (mg C m day ) from 1999 to 2009, (b) simulated growing season (day of year 120 to 330) daily CH 4 flux (mg C m -2 day -1 ) from 2004 to 2008 (blue area indicates the variation range as recoreded for the Mer Bleue Bog averaged from 12 collars (Moore et al. 2011, Fig. 2), (c) relationship between the simulated growing season daily and (d) annual average CH 4 flux (mg C m -2 day -1) in log 10 scale and the simulated daily water table depth (m) from 1999 to 2009.

10FigureFig. 10 .
Figure 10.Carbon budget (a) and nitrogen budget (b) for the Mer Bleue peatland based on simulated averages from 1999 to 2004.Pools are in g m -2 and fluxes are in g m -2 year -1 .
(+) indicates a positive relation between the change in the parameter and the change C and N pools or fluxes.(-) indicates an inverse relation between the change in parameter and the change in C and N pools or fluxes.

Figure 11 .Fig. 11 .
Figure 11.Simulated C pools (a), N pools (b) and NPP (e) in plants and PFTs, N uptake rate by in plants and sequestration rate in peat (c), N absorbed by moss and N content in mosses (d), N mineralization rate and N output from peat (f) at annual wet N deposition of 1.5g N m -2 yr -1 from 1999 to 2039.

Table 1 .
State variables in the model.
* Units were standardized to 1 m 2 area of Peatlands for model output.

Table 4 .
Assumed and calibrated parameters.
, the average annual GPP was 555 gC m −2 yr −1 , of which 70 % was contributed by shrubs and 26 % by mosses.Average annual ER was 526 gC m −2 yr −1 , 73 % of which was emitted from the soil surface produced in HR of microorganisms and AR in roots.The difference of GPP and ER resulted in 286 gC m −2 yr −1 of NPP of plants on average, whereas the www.geosci-model-dev.net/6/1173/2013/Geosci.Model Dev., 6, 1173-1207, 2013

Table 5 .
Results of sensitivity analyses.The values shown are the average relative changes in model output per change of parameter

Table 6 .
Observed (Obs.), simulated (Sim.), and the difference (D) between observed and simulated annual GPP, ER and NEE (units: gC m −2 yr −1 ) and summer average water table depth (unit: m) from 1 May to 30 October for 6 yr for the Mer Bleue Peatland.