CO 2 ﬂuxes and ecosystem dynamics at ﬁve European treeless peatlands – merging data and process oriented modeling

. The carbon dioxide (CO 2 ) exchange of ﬁve different peatland systems across Europe with a wide gradient in land use intensity, water table depth, soil fertility and climate was simulated with the process oriented CoupModel. The aim of the study was to ﬁnd out whether CO 2 ﬂuxes, measured at different sites, can be explained by common processes and parameters or to what extend a site speciﬁc conﬁguration is needed. The model was calibrated to ﬁt measured CO 2 ﬂuxes, soil temperature, snow depth and leaf area index (LAI) and resulting differences in model parameters were analyzed. Finding site independent model parameters would mean that differences in the measured activity, mineralization–immobilization, nitriﬁcation, denitriﬁcation and plant respiration t Q10bas base temperature for the microbial activity, mineralization–immobilization, nitriﬁcation and denitriﬁcation at which the response is 1


Introduction
In recent years, many data sets have been collected from a number of sites and across multiple years, containing detailed and high-resolution measurements of carbon (C) fluxes, plant and soil characteristics, meteorological and water table data (Baldocchi et al., 2001;Baldocchi, 2007).
Several of the measured sites are peatlands, which have accumulated a vast amount of C since the last deglaciation. Under drained conditions, peatlands have a high carbon dioxide (CO 2 ) emission potential (e.g., van den Bos, 2003;Lohila, 2004;Drösler et al., 2008;Maljanen et al., 2010). Understanding the processes driving CO 2 emissions is essential in the development of management practices to reduce greenhouse gas emissions.
Direct comparison of measured data can be used to explore the effect of single variables if the site conditions are similar or differ only in few variables, e.g., in manipulation experiments (Chivers et al., 2009;Ward et al., 2013) or different vegetation types at the same site (e.g., Chojnicki et al., 2010). However, the sites in this study have very different characteristics with respect to climate, hydrology, current and former land management, vegetation and soils. Direct site comparisons of measured flux data (e.g., Alm et al., 1999;Humphreys et al., 2006;Lund et al., 2009;Drewer et al., 2010) are often uninformative when trying to distinguish between responses of several individual factors. Typically, multiple factors are linked and interact with each other complicating the analysis. Therefore, important drivers at one site might not play a significant role on another site (e.g., Lafleur et al., 2005). Process oriented modeling provides a method to identify to what extent observations at different sites can be described by the same processes, while accounting for such interactions.
Process oriented modeling requires (1) that the model can describe the observations and (2) that the parameters used in the model to describe the observations can be estimated from available data. Typically, studies focus on demonstrating how well the model can describe a certain set of data (e.g., van Huissteden et al., 2009;Calanca et al., 2007;Frolking et al., 2001;St-Hilaire et al., 2010). In contrast, the focus of this study was exploring differences between the sites while model performance was subordinate. Process oriented models often require a large number of input parameters which are usually difficult to estimate based on available data from less intensively investigated sites (Juston et al., 2010). Parameters may interact with each other and the available information does not allow for a single or unambiguous mathematical solution (Beven and Freer;2001, Beven, 2006van Oijen et al., 2013). However, for all sites in this study, accurate gas flux measurements in combination with detailed measurements of soil and plant conditions were available. Such extensive measurements have been demonstrated to be useful in identifying the governing properties for specific sites. For example, the modeling of CO 2 from forest sites has shown that dynamics of CO 2 fluxes are restricted to a certain range of parameter values (Wu and Jansson, 2013;Wu et al., 2013a).
A systematic evaluation of one model against data from multiple sites with a common set of parameters will allow for a better understanding of processes not only at the individual sites but also on the site specific differences which control the resulting fluxes (e.g., Calanca et al., 2007;van Huissteden et al., 2006;van Huissteden et al., 2009). This is a necessary precondition for accurate predictions of CO 2 fluxes under different climate scenarios or at different locations. On peatlands, some attempts have been made to consider site differences using simplified process models on national (e.g., ECOSSE; Bell et al., 2012) and global scales (e.g., InTec; Ju and Chen, 2005;McGill;St-Hilaire et al., 2010) and up to a millennial timescale (Schuldt et al., 2013). However, we are not aware of any studies comparing differences in parameter distributions of CO 2 -related processes between treeless peatland sites using an uncertainty-based approach and a detailed process oriented model running on site scale.
In this work the CoupModel was used, which is a detailed process oriented model coupling heat and mass transfer for soil-plant-atmosphere systems (Jansson and Karlberg, 2010). The CoupModel was chosen for the following reasons: the model was designed for a wide range of soil types and different ecosystems and applications (see  for review) which might be useful as some of the sites in this study are already quite degraded and might not respond like a typical, intact peatland anymore. The model has been shown to be capable of simulating all three main greenhouse gases from peatlands: CO 2 , nitrous oxide (N 2 O) (Norman et al., 2008) and methane (CH 4 ) (Ravina, 2007). Furthermore, the CoupModel includes detailed submodules for the most relevant processes in the carbon cycle: it predicts plant growth, plant transpiration and autotrophic respiration, soil nitrogen (N) and C processes, energy and heat fluxes, soil temperature, soil frost and snow depth. It supports an hourly time step for input and output data and can run in even finer time resolution, which is necessary for analyzing e.g., chamber flux data. The user can select between different submodels, different equations and different complexities and easily access all parameters via a user interface. Calibration procedures with randomized parameter values and methods for visualization and detailed analysis of the model output are supported. An extensive model description can be found in Jansson and Karlberg (2010). The model and its documentation as well as several tutorials for its application can be downloaded from the CoupModel home page (CoupModel, 2014). The main aim of this study was to find out to what extend the large differences in measured CO 2 fluxes between five data rich European flux measurement sites can be solely explained by the differences in meteorology, water table and management. Therefore, the process oriented CoupModel was applied using an uncertainty-based Monte Carlo approach. Specific objectives were I. to identify differences and similarities between various sites in CO 2 -related processes, corresponding parameters and responses to forcing data; II. to identify and discuss the impact of available data for estimating key parameters in CO 2 flux models in general; III. to identify problems related to the model representation of the different ecosystem processes for open peatlands.

Description of sites and investigations
The CoupModel was applied to five treeless peatland sites with a wide gradient in land use intensity, water level, soil nutrient status and mean annual temperature (Table 1). Together with the climatic gradient from northern Finland to southern Germany and a different growing season, this leads to great differences in amplitude and dynamics of gross primary productivity (GPP), ecosystem respiration (R eco ) and different amounts of biomass. This is reflected in the annual accumulated net ecosystem exchange (NEE) based on measurements, ranging from −395 to 636 g C m −2 (Fig. 1).
Dynamic forcing data for model input (water table and meteorology) were available from measurements at all sites (Table S1). Data used for model parameter constraint included measurements of LAI, soil temperature and NEE (Table S2). Measured NEE was partitioned into R eco and GPP by the use of empirical models. At sites where the eddy covariance method was used, R eco was derived from nighttime NEE, otherwise it was taken from opaque chamber measurements.. The empirical R eco models are based on temperature (Lloyd and Taylor, 1994), while light-level-based functions were used for GPP according to Falge et al. (2001). Corrections and gap filling at flux tower sites was done according to the methods described in Reichstein et al. (2005). A detailed description is given in the references listed in Table S2. Though R eco and GPP are not explicitly measured, this will be called measured data in the following for simple distinction from the simulated fluxes by the CoupModel.
The northernmost site, Lompolojänkkä fen (Lom), located in Finland is a nutrient-rich natural mire with sedges, shrubs and mosses. Mean air temperature from 2006 to 2010 was −1.4 • C and the mean groundwater table during the snowfree season was close to the peat surface. Data for model calibration were available from 2006 to 2010 and consisted of eddy covariance (EC) and automatic chamber data of CO 2 fluxes, snow depth and leaf area index (LAI) measurements. A detailed description of the site and measurement methods can be found in Aurela et al. (2009), Drewer et al. (2010 and Lohila et al. (2010).
The Scottish site, Auchencorth Moss (Amo) is an ombrotrophic bog, with vegetation consisting of grasses, sedges and soft rushes, covering a primarily Sphagnum base layer. The site is managed for low-intensity sheep grazing with less than one livestock unit per hectare, but this was not accounted for in the model. Amo encompasses a small area of peat extraction in the southwest of the catchment, which is unlikely to fall within the flux footprint of the EC system. The site was drained over a century ago, however, the drains are no longer considered to be in operation. The mean water table was −12.5 cm between 2006 and 2010. Mean temperature during this period was 10 • C, CO 2 data from EC during the same period was used for model calibration. A detailed description of the site and measurements can be found in Helfter et al. (2014), Drewer et al. (2010) and Dinsmore et al. (2010).
Horstermeer fen (Hor) is located in the Netherlands in a drained natural lake. It used to be agricultural land but was abandoned more than 15 years ago. The water table was raised during restoration leading to a mean value of −10 cm during the simulation period from 2004 to 2010. It became a seminatural grassland, a nature reserve without any mowing management. The vegetation is very heterogeneous with reeds, grasses and small shrubs . The mean temperature during the simulation period was 10 • C. CO 2 fluxes were measured half-hourly by EC and biweekly with opaque chambers between 2004 and 2010. A detailed description of the site and measurement methods can be found in Hendriks et al. (2007).
Freisinger Moos (FsA and FsB) is a drained nutrient-rich fen in the south of Germany. The two sites FsA and FsB lie next to each other in a drained sedge meadow which was cut once per year. The mean annual hay yield was 4.19 or 4.07 t dry weight ha −1 a −1 for FsA and 5.67 or 6.17 t dry weight ha −1 a −1 for FsB for the years 2010 and 2011, respectively. FsB is located in a small depression with a mean water level of −20 cm compared to −25 cm for FsA during The measured R eco and empirical modeled GPP during the measurement period of each measurement day were used for parameter constraint, empirically modeled values between measurement days were only used for visualization and comparison.

Model description
CoupModel v4 from 12 April 2013 was used for simulations. The current version can be downloaded from the CoupModel home page (CoupModel, 2014). A detailed description can be found in Jansson and Karlberg (2010). The model represents the ecosystem by a description of C and N fluxes in the soil and in the plant. It includes all main abiotic fluxes, such as soil heat and water fluxes, that represent the major drivers for regulation of the biological components of the ecosystem. The most important equations with the corresponding parameters and switches differing from the default setup in the used version can be found in Tables S3, S4, S5 and S6. The major model assumptions relating to the model application to peatlands are described below. Figure 2 shows a scheme of the main carbon fluxes and pools in the current CoupModel setup. gross primary production k gresp growth respiration coefficient k h rate coefficient for the decay of the slow C pool k l rate coefficient for the decay of the fast C pool k tot total rate of decomposition calculated from k h , k l and SOC of the corresponding pools in the upper 30 cm k mrespleaf maintenance respiration coefficient for leaves k rn extinction coefficient in the Beer law used to calculate the partitioning of net radiation between plant canopy and soil surface LAI leaf area index ME mean error m retain coefficient for determining allocation to mobile internal storage pool N nitrogen NSE Nash-Sutcliffe efficiency p ck speed at which the maximum surface cover of plants is reached p θp power coefficient in the response function of microbial activity in dependency of soil moisture p θSatact activity under saturated conditions in the soil moisture response function for microbial activity, mineralizationimmobilization, nitrification and denitrification p θUpp water content interval in the soil moisture response function for microbial activity, mineralizationimmobilization, nitrification and denitrification R 2 coefficient of determination R eco ecosystem respiration SOC soil organic carbon T amean assumed value of mean air temperature for the lower boundary condition for heat conduction T aamp assumed value of the amplitude of the sine curve , representing the lower boundary condition for heat conduction T MatureSum temperature sum beginning from grain filling stage for plant reaching maturity stage T DormTh critical air temperature that must be undershot for temperature sum calculation T EmergeSum air temperature sum that is the threshold for start of plant development T EmergeTh critical air temperature that must be exceeded for temperature sum calculation t Q10 response to a 10 • C soil temperature change on the microbial activity, mineralization-immobilization, nitrification, denitrification and plant respiration t Q10bas base temperature for the microbial activity, mineralization-immobilization, nitrification and denitrification at which the response is 1 ε L radiation-use efficiency

Meteorological driving variables and integration time step of the model
Hourly values of global radiation, relative humidity, precipitation, wind speed, and air temperature, measured at each site were used as input. Data was gap filled by simple linear interpolation for gaps < 6 h. Larger gaps were filled by values from other adjacent climate stations. At Hor the station used for gap filling provided only daily values. Hourly values were retrieved assuming uniform distribution over 24 h for precipitation, wind speed and relative humidity and sinusoidal distribution for temperature and global radiation. Model performance was only evaluated for the years when meteorological data was available. The simulations were started 2 years prior to the evaluation period, so the system (in particular the plant) could adapt to the site conditions and become more independent of initial values. Data from the available years was copied to previous years if not available from an adjacent climate station.
The model's internal time step was half-hourly for abiotic processes and hourly for nitrogen-and carbon-related processes.

Dynamic coupled heat and water model for above soil surface conditions
An interception model for both, radiation and precipitation, a snow model and a surface pool model was used to provide boundary conditions at the soil surface. Interception and plant evaporation was dependent on the simulated leaf area index of the plant as well as the degree of coverage, while transpiration depended additionally on the simulated water uptake of the plant. Cloud fraction was calculated from global radiation input and latitude. Incoming radiation was partitioned between one part which was absorbed by the plant canopy and another part which reached the soil. Surface temperature was simulated based on an energy balance approach, where the radiation reaching the soil equals the sum of sensible and latent heat flux to the air and heat flux to the soil. Soil evaporation was derived from an iterative solution of the soil surface energy balance of the soil surface, using an empirical parameter for estimating the vapor pressure and temperature at the soil surface. Vapor pressure deficit was calculated from the relative humidity input. Snowfall was simulated from precipitation and air temperature, snowmelt from global radiation, air temperature and simulated soil heat flux. Surface runoff was controlled by a surface pool of water that covers various fractions of the soil surface. Under oversaturated periods the flow of water in the upper soil compartment could be directed upwards, towards the surface pool. Surface runoff was calculated as a function of the amount of water in the surface pool.

Dynamic heat and water model for the soil
The soil profiles were divided into 12 layers with an increasing layer depth from 5 cm for the upper layer to 100 cm in the lowest layer. Heat flow between adjacent soil layers were calculated based on thermal conductivity functions accounting for the content of ice and water. The heat flow equation is based on a coupled equation accounting for the freezing and thawing in the soil (Jansson and Halldin, 1979). Convection was not accounted for. The lower boundary was calculated as temperature based on a sine variation at the soil surface and a damping depth for the whole soil profile as well as a parameter for the annual mean temperature T amean and annual amplitude of temperature T aamp at the site (a list of symbols and abbreviations can be found in Table 2). Soil water depended on infiltration to the soil, soil evaporation, water uptake by plant roots and groundwater flow. Soil moisture, represented as liquid water content, was calculated based on the water storage and temperature. Water flows between adjacent soil layers were calculated according to the Richards equation (1931), considering hydraulic conductivity, water potential gradient and vapor diffusion. Soil water characteristics were described by the Brooks and Corey (1964) equation between two threshold water tensions, while a log-linear expression was applied at higher water tensions and a linear expression at water contents close to saturation. Unsaturated conductivity was simulated according to Mualem (1976), additionally accounting for the conductivity in macro pores. The groundwater level was defined by assuming a continuous zone of saturation from the water table level down to the lower boundary of the considered soil pro-file. To force saturation at the measured groundwater level, water was added to or removed from the corresponding layer.

Vegetation
Vegetation was simulated according to the explicit big leaves concept (e.g., Dai et al., 2004) but only one plant canopy layer, representing the complete plant community, was defined. Albedo, LAI, vegetation height and vegetation cover were simulated. Permanent, perennial vegetation was configured with maximal plant height of 0.6 m, a lowest root depth of −0.6 m and a maximal plant cover of 100 %. Grain development was assumed to play a minor role and was therefore disabled. Plant respiration was assumed to be dependent on growth and maintenance (e.g., Hansen and Jensen, 1977).
For leaf assimilation, the light-use efficiency approach (Monteith, 1972;Monteith and Moss, 1977, see e.g., Hilker et al., 2008 for review) was used, at which total plant growth is proportional to the global radiation absorbed by canopy but limited by unfavorable temperature and limited soil water. For simplicity, plant assimilation was simulated independent of dynamics in N availability. This might be justified as none of the sites were fertilized in the recent years and the vegetation community was assumed to be adapted to the nutrient conditions at each site. Differences in N availability between sites are included in the radiation efficiency (ε L ). Plants were assumed to be well adapted to wet conditions (Keddy, 1992, Steed et al., 2002, including aerenchyma to tolerate water saturated soil conditions (Jackson and Armstrong, 1999). Plant stress due to high water saturation was therefore disabled.
Plant development started every spring when the accumulated sum of air temperatures above a threshold value (T EmergeTth ) reached the value of T EmergeSum . Both parameters were calibrated (Table S4). The accumulation of temperatures started when the day length exceeded 10 h. Snow cover hindered shooting by reducing the radiation passing through to the plant, while low soil temperatures reduced plant water uptake.
Besides a small amount of litter fall occurring during the whole plant growth period (Robson, 1973;Duru and Ducrocq, 2000;Fulkerson and Donaghy, 2001), senescence was assumed to start after the plant reached maturity and therefore depended on growth stage (e.g., Thomas and Stoddart, 1980) and temperature sums (e.g., Davidson and Campbell, 1983). As this was not yet directly supported by the model, the stem pool was used for brown, senescent, standing biomass. Therefore new assimilates were constantly allocated to roots and leafs only, while existing leaf biomass was reallocated after maturity to the stem pool. A third stage of litter fall was configured depending on a minimum threshold temperature sum: five consecutive days in the autumn with day lengths shorter than 10 h and with temperatures below T DormTth terminated the growing season and plants went into dormancy.
During litter fall part of the C is stored in a mobile pool, which can be then reused for regrowth in spring in the next year (e.g., White, 1973;Wingler, 2005).
Harvest took place at FsA and FsB. Based on observations in the field, 85 % of the aboveground plant material was removed at harvest. Harvest dates were known and implemented in the model. After harvest the growth stage was allowed to be reset to a lower value (e.g., Thomas and Stoddart, 1980). Reallocation of C from roots to leaves could take place as reported for, e.g., Festuca pratensis (Johansson, 1992;1993).

Soil carbon and nitrogen
The organic substrate was represented by two C and N pools for each of the 12 soil layers: one with a slow and one with a high turnover rate coefficient. Decomposition products from the fast pools are partitioned into CO 2 which is released to the atmosphere and C which is partly moved to the slow pools and partly returned to the fast pools. Decomposition products from the slow pools are partly released as CO 2 and partly returned to the slow pools. The initial values for the amount of C and N per layer was given by measurements and partitioned into the two pools for each layer according to the measured C : N ratio as described in Sect. 2.2.5 and Table 3. Beside the turnover rate coefficients and amount of substrate in each pool per layer, decomposition rates depended on the response to soil moisture and temperature in the corresponding layer.
As the rate coefficients for decomposition were expected to strongly affect each other, only the coefficient for the fast decomposition pools were calibrated. The coefficient for the slow pools (k h ) was kept constant at a low value of 2 × 10 −8 day −1 during the calibration runs, which might be justified as decomposition of resistant carbon is less responsible for the variation in soil respiration (e.g., Whalen et al., 2000).
Nitrogen-and methane-related processes were considered by a model including the most important pathways and fluxes. However no emphasize on the calibration of these processes was made in this study since the current objective was on CO 2 fluxes from the peatlands.

Independent approach to find values of site specific parameters
Dry and wet N deposition, latitude and thickness of the organic layer were used as constant site specific input. Water retention parameters were assigned to each soil layer according to soil data from each site. However, at Amo and Lom, water retention and, at all sites, unsaturated conductivity was assigned from the CoupModel soil database as suggested by Lundmark (2008) for peat soils. Measured total soil organic carbon (SOC) per layer was partitioned to the two SOC pools per layer on the basis of the measured total C : N ratio per layer whereas the initial C : N ratios of the slow decomposing pools were assumed to be 10, while for the fast pool's 27.5 was chosen according to measured C : N of leaf tissues at FsA and FsB (Table 3).

Parameter calibration approach
The aim of the calibration was to find out to what extent the same parameter values could be used for all sites compared to a site specific representation. A stepwise approach was carried out starting with finding the best site specific parameter representations and then trying to merge them to common values valid for all sites. Finally, the common representation was revised to some few parameters showing great site specific effects on model performance. An overview of the different steps can be found in Fig. 3, details on the calibration procedures are presented as supplementary material.
For the basic calibration (step I, Fig. 3), 350 000-700 000 runs were performed for each site. 45 parameters, which were suspicious of eventually being site specific, were selected and calibrated with an assumed uniform random range  ( Table S4). Parameter ranges were then constrained based on selected runs (steps I and II; Fig. 3), showing acceptable performance to multiple variables (Table S7), measured at the sites. Several additional multiple calibration runs were performed, with few selected parameters each, to unravel parameter interactions (step III; Fig. 3). A number of simulations were also made by single value representations of parameters (step IV; Fig. 3) to visualize the impact of certain parameter values on interacting parameters and on performance. These runs are called single runs in the following, numbered C1-C7 and described in Table S8.
Selection of runs and evaluation of performance was based on three indices: coefficient of determination (R 2 ) assesses how well the dynamics in the measurement derived values are represented by the model. Mean error (ME), also called y intercept (Willmott, 1982) indicates a lag or lead between model predictions and measured data (Moriasi et al., 2007). Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) accounts for both deviation of dynamics and magnitude. It ranges from −∞ to 1, whereas 1 means the best fit of modeled to measured data and values < 0 indicates that the mean measured value is a better predictor than the simulated value, which indicates an unacceptable performance (Moriasi et al., 2007).

Model performance -results of basic calibration and selected common configuration
Model performance showed distinct differences between the sites, depending on the investigated variable and on the num-ber of considered runs (Table 4). Figure 4 shows the differences between measurements and model C1.

Fluxes
At all sites the dynamics in R eco fluxes were simulated considerably better than GPP (Table 4). Performances for NEE were worse as simulation errors in GPP and R eco are summed up.
With respect to R eco and GPP, the selected single runs represent a parameter configuration close to the best ones possible in the tested range: their R 2 value did not differ by more than 0.05 from the best result achieved in the multiple calibration, while ME values were smaller than 0.1 g C m −2 day −1 . Clearly, the lower R 2 and higher ME values in single runs for biomass and LAI simulation indicate that none of the runs could give the best results for all variables at the same time. For example, the best values for GPP can only be achieved if poorer performance would have been accepted for other parameters such as winter R eco or LAI (see criteria for accepted runs in Table S7).
The ME values in Table 4 show a clear overestimation of winter fluxes by 3.21 and 2.11 g C m −2 day −1 for the single runs at FsA and FsB, respectively, and a weaker overestimation for the accepted runs. The overestimation was less pronounced at Amo (0.13 g C m −2 day −1 ) and Lom (0.01 g C m −2 day −1 ). At Hor, winter fluxes were underestimated with a ME of −0.26 g C m −2 day −1 . This was reflected in the accumulated NEE (Fig. 4) leading to a much higher CO 2 loss compared to the CO 2 balance estimated by the empirical model approach at FsA and FsB. At Lom, higher accumulated NEE due to the overestimation of winter R eco was visible in the first months of each year. It was nearly compensated due to the underestimated spring R eco or overcompensated due to GPP overestimation as, e.g., in summer 2006, which was very dry.

Explanatory variables
Of all variables, the highest R 2 values were achieved for soil temperature at all sites. Temperatures in deeper soil layers (−50 or −60 cm) had better fits than in upper layers with R 2 values close to 0.9 or higher and a maximum mean deviation of 0.15 • C. The fit of modeled vs. measured snow depth, which was only available at Lom, had a R 2 value of 0.75 with a mean error of less than 10 cm. Simulation of LAI represented the measurements quite well with R 2 values of between 0.53 and 0.76 and a mean error of maximum 0.12 m 2 m −2 . An exception was Hor, where LAI was underestimated by a ME of −0.61 and −1.49 m 2 m −2 in the accepted 75 runs and in the selected single run C1, respectively. At Hor, root biomass was underestimated in the single run by a ME of −281 g C m −2 and living leaf biomass by −122 g C m −2 .
In most of the runs of the basic calibration at Hor, either GPP was overestimated or leaf biomass and LAI were underestimated. Therefore, beside the common configuration C1, a different configuration was tested where plant respiration and litter fall parameters for Hor were set to much lower values than in the tested range to fit to GPP and LAI at the same time. However, this reduced performance for R eco R 2 to 0.66 compared to 0.75 in C1 and led to an overestimation of winter R eco with a ME of 0.75 g C m −2 day −1 .

Parameter constraint
Site specific calibration was needed for the speed at which the maximum surface cover is reached (p ck ), the mean value in the analytical air temperature function (T amean ), temperature sum for reaching plant maturity (T MatureSum ), coefficient for determining allocation to mobile internal storage pool (m retain ), decomposition rate of the fast SOC pools (k l ) and radiation-use efficiency (ε L ).
Activity under saturated conditions (p θ Satact ), threshold temperature for plant dormancy (T DormTh ), response to a 10 • C soil temperature change on the microbial activity (t Q10 ) and base temperature for the microbial activity (t Q10bas ) covary with performance indices but showed different patterns for different validation variables and for different sites.
Most of the parameters did not show any influence on performance indices within the tested range (Fig. S1), demonstrating that either the relatively low effect of the parameter was overcompensated by the effect of more sensitive parameters, or the range used for calibration is sufficiently constraining. Each of these parameters did not reduce model performance indicated by R 2 by more than 0.05 for GPP or R eco after setting them to a common value.

Correlations between parameters
In the basic calibration, the following parameters were identified to interact with other parameters: p ck covaried with the extinction coefficient in the Beer law (k rn ) which is used to calculate the partitioning of net radiation between canopy and soil surface. Strong, linear, negative correlation between coefficients for growth (k gresp ) and maintenance respiration (k mrespleaf ) was detected.
The effect of the different parameters in the water response function, p θ Satact , p θUpp and p θp , compensated each other. They could not be constrained without a very high measurement resolution of fluxes and water table combined with high water table fluctuation at the same time. Therefore, p θ Upp and p θ p were set to default values and p θSatact was con-strained by additional multiple runs together with k l . Differences between sites in k l are reduced with higher p θ Satact values (Fig. 5), however, higher p θSatact values increase overestimation of winter R eco at FsA and FsB (Figs. 6, 7d). A wider range of p θ Satact values was acceptable for summer R eco (Fig. 6).
Besides moisture response, decomposition rate (k l ) and temperature response (t Q10 , t Q10bas ) control soil respiration. The effect on R eco was cofounded by plant respiration. Different patterns for different sites and variables for each of the parameters were even more pronounced when only k l , t Q10 and k mrespleaf were in calibration (Fig. 6).
Single runs with different configurations (Fig. 7) revealed that higher plant respiration as well as steeper temperature response can lead to less overestimation of respiration in winter Figure 5. Dependencies between the parameters for decomposition rate and saturation activity for the different sites, based on additional multiple runs.
( Fig. 7d) but lead to reduced performance (Fig. 7c). In all single runs, despite the different configurations, FsA always showed the highest k l , while Amo had the lowest (Fig. 7a). A higher saturation activity reduces the difference in k l values, but leads to higher overestimation of winter fluxes.

Discussion
Despite not being specifically developed for peatlands, the CoupModel was able to simulate measured fluxes quite well. The model was run in a simple configuration with only two SOC pools per layer, no explicit representation of microbes, and only one plant layer. Even though the CoupModel was capable of adequately reproducing the measurements. Several points were identified where further peatland specific processes or more detailed representations might improve the model. Those are discussed in the following subsections.
From the 45 calibrated parameters, 8 parameters could be identified to actually need a site specific representation to achieve acceptable performance. Those parameters are discussed in Sects. 4.3-4.10. The remaining 37 parameters were not sensitive in the tested ranges, even though site specific values could have been expected; for example, it is known that grassland species differ in their assimilation and growth response to temperatures (Billings et al., 1978;Wohlfahrt et al., 1999). Plant respiration rates in graminoids differ between species (Poorter et al., 1991;Scheurwater et al., 1998;van der Werf et al., 1988) and depend, among other factors, on light (Rovira, 1969;Bahn et al., 2013), nutrient (Paterson and Sim, 2000) and moisture conditions (Crow and Wieder, 2005) as well as on cutting regime (Bokhari, 1977). Allocation fractions to different plant parts differ between species and depend on nutrient conditions (Aerts et al., 1991;Gong et al., 2014) as well as shading (Bahn et al., 2013). Values for specific leaf areas are species specific (e.g., Poorter et al., 1990;Reich et al., 1998) and differ in response to nutrient availability (Meziane and Shipley, 1999). Leaf lifetime (e.g., Ryser and Urbas, 2000) as well as leaf and root turnover rates (Schläpfer and Ryer, 1996) vary between graminoid species.
The five peatland sites largely differed in their vegetation composition, plant life-forms and species. Nevertheless, common values for all sites could be applied for parameters related to these processes, without reducing model performance on R eco and GPP in R 2 values by more than 0.05. That shows that either the studied sites on a vegetation community level did not differ much in these processes, or that the impact of those parameters is subordinate compared to the impact of other parameters, meteorological input and other site conditions. Therefore, models with a focus on multiple-year carbon fluxes do not need a site specific interpretation of these processes.

Model initialization
Many models use spin-up routines of many years until SOC pools reach a steady state (e.g., Dimitrov et al., 2010;Smith et al., 2010;Thornton and Rosenbloom, 2005). Here, measured C : N values were used to partition the SOC between pools, while ranges for parameter values were chosen in a way that the amount of carbon in the soil pools did not change very drastically. However, no further effort was made to force the pools to be in equilibrium. It was assumed that this might not be the case in the real world either: drainage ditches at FsA and FsB are still maintained, leading to high carbon losses and changes in substrate quality. Land use at Hor was quite recently changed from a fertilized and deeply drained cropland to a nature reserve with a restored water table. Additionally, Amo used to be more intensively managed and drained, but the drainage system was not maintained. Land use history was not known and SOC measurements were available from only one date per site. The measured carbon fluxes were therefore the only indication about carbon loss or addition to the complete system, while changes in relative pool sizes were not known. The partitioning of the SOC has implications on the parameter distribution for the rate coefficient for decomposition, which is discussed in Sect. 4.10.

Model performance
The best performance achieved highly differed between the different validation variables and between the different sites. This was not only caused by the model's ability to simulate the different output parameters but also due to measurement quality, measurement uncertainty, measurement methods (temporal and spatial resolution) and heterogeneity of the sites.   Table S8).
GPP was simulated markedly poorer as compared to R eco at all sites and not only in the single runs, but also in the complete set of performed multiple runs. An explanation might be that in the model the whole plant community consisting of different individuals, species and even functional types, with different life cycles and adaptations to light availability and temperature was simplified to only one plant. Especially mosses, which differ largely from vascular plants in respect to their ecology and response to water, temperature and light conditions (Gaberščik and Martinčič, 1987;Harley et al., 1989;Murray et al., 1993;Turetsky, 2003), might be important at the moss-rich Lom and Amo sites. The vegetation at Hor consists of species with very different strategies and requirements for nutrients and water. At FsB, reeds, which are known for a late emergence, were well present in some of the years and hardly appeared in other years. FsA is relatively species-rich and several of these species are abundant only during parts of the vegetation period. Also, using a more complex photosynthesis model such as in e.g., Farquhar et al. (1980Farquhar et al. ( , 2001 and testing a wider range of parameters might lead to a better fit. Including plant stress due to high water levels and nutrient limitation might improve the performance on some sites. For example, Sagerfors et al. (2008) found photosynthesis to be limited also by too high water levels, so that the McGill wetland model assumed reduced photosynthesis outside a water level range of −10 to −20 cm (Wu et al., 2013b). Furthermore, GPP cannot be measured directly neither by the chamber nor the EC method. Instead it was derived from NEE and R eco or nighttime NEE, including the uncertainty of two different measurements and empirical modeling.
Heterogeneity of vegetation was very distinct at Hor, which might explain the difficulties to simulate the right amounts of GPP and biomass at the same time. The biomass and LAI taken into account for this study might not be fully representative of the whole EC fetch for all wind directions. Hor is also a site which deviates strongly with respect to other sites, with recent large changes in management. It is in successional transition from intensively used dairy farming meadow (approximately 20 years ago) towards reed fen with willow thickets. Soil and vegetation still show the imprint of high nutrient levels derived from manuring practices (e.g., patches with abundant Urtica dioca). This likely still affects GPP. These features could be a better explanation of the deviating GPP than the additionally tested configuration with strongly reduced litter fall and plant respiration rates.
Even though the winter fluxes are small compared to the summer fluxes, they have a marked role in the annual NEE balances (Fig. 4). Overestimation of winter R eco in combination with slightly underestimated winter GPP leads to high overestimation of annual accumulated NEE, emphasizing the importance of winter flux dynamics in the annual balances. At all sites except Hor, winter R eco was overestimated in the selected single run. For FsB and especially FsA this was also true for all multiple runs. As R eco at Lom and Amo are typically relative low, the effect was less pronounced.
Several different reasons for the winter R eco overestimation are possible: explanations due to model setup and parameterization are discussed in the Sects. 4.7, 4.8 and 4.9. Additionally, gases might be trapped within the snow and under the ice Maljanen et al., 2010) and therefore captured by the measurement instruments only in springtime, when they are released. A gastight ice cover was not realized in the current model setup. Frozen or ice covered soils are quite common at the boreal Lom, but also at FsA and FsB which have a more continental climate than the other sites.
The ability of the model to simulate soil moisture could not be evaluated, as this variable was measured only at Lom, where the soil was close to saturation throughout the year. Therefore, and as groundwater level was used as input, hydraulic properties could not be constrained. Furthermore, swelling, shrinking and hysteresis effects which are important factors in hydraulic characteristics of peat soils (e.g., Kellner and Halldin, 2002) were not accounted for. This could have an effect on model performance and parameter values, especially those related to the soil moisture response.

Soil temperature dynamics
Due to the isolating impact of the snow cover (e.g., Zhang, 2005), the value of mean annual soil temperature (T amean ) was expected to be slightly higher than the mean annual air temperature. Constrained values of soil temperature were 1.5-5 • C higher than the mean annual air temperature at all sites. If the model was run under different conditions without further fitting, factors causing differences between mean annual soil temperatures and corresponding air temperature need to be considered.

The role of soil temperature and GPP to constrain the plant cover
Accepted fits for soil temperature in the uppermost measured soil layer led to p ck values close to the measured coverage of vascular plants for each site. Therefore, the measured coverage could directly be used in the configuration C1 (Fig. 6a).
Setting p ck to a common value of 100 % reduced the differences in ε L between the site's C7 (Fig. 6e), but led to underestimation of soil temperature in the uppermost soil layer by at most −0.45 • C in ME at Amo. An explanation could be that mosses are contributing to the plant coverage in respect to GPP but not to temperature, especially at sites where they are the main peat-forming material.

Start of senescence
Site specific calibration was needed for the temperature sum initiating the start of senescence (T MatureSum ). However, if the resulting day of the year was plotted instead, the differences between sites became small (Fig. 6) and setting it to the mean value of all sites did not reduce model performance in GPP R 2 by more than 0.05. Induction of senescence with graminoids is known to depend on both temperature and day length (Nuttonson, 1958;Proebsting et al., 1976;Thomas and Stoddart, 1980;Davidson and Campbell, 1983). However, the differences between the sites in this study could be explained solely by the relative day length. The proportion of C in the plant which does not become litter, but instead is stored for shooting in the next year (m retain ), differed largely between sites. At Lom, a value of at least 40 % led to and acceptable performance while a maximum of 3 % was found for FsA and FsB; a mean value of 20 % would reduce R 2 of GPP by at least 0.04 for these sites. At Amo and Hor neither a value of 3 % nor 40 % reduced R 2 of GPP by more than 0.01. An explanation for the low m retain at FsA and FsB could be that the same pool is used for regrowth after being cut and therefore not available for shooting anymore, as the regrowth rate in both early spring and after being cut depend on the carbohydrate reserve (White, 1973;Davies, 1988;Klimeš and Klimešová, 2002). Steele et al. (1984) conclude that defoliation late in the year will affect spring regrowth. At Lom, high m retain might be an adaption to the short vegetation period (Kistritz et al., 1983). Evergreen parts of the vegetation like dwarf shrubs, lower leaf parts of graminoids and mosses were not accounted for, which also affects regrowth in spring. Saarinen (1998) found that 60-70 % of shoots and 20 % of green biomass in a Carex rostrata fen survived the winter and, therefore, hypothesized based on comparisons with other studies that the proportion increases with increasing latitude.
The storage pool is an important parameter needing site specific calibration but can be fitted if several measurements during spring and early summer of either GPP, biomass or LAI are available.

Radiation-use efficiency
As plants were not nutrient limited in the model setup, the lowest values for ε L were expected under the most nutrientpoor conditions (Longstreth and Nobel, 1980;Reich et al., 1994;Haxeltine and Prentice, 1996;Gamon et al., 1997;Wohlfahrt et al., 1999). The opposite was true if site specific values were used for p ck . However, a common value for p ck reduced the differences in ε L and led to low ε L at the ombrotrophic Amo site, but to an even lower value at the minerotrophic Lom site. The nutrient status of the soil can therefore not explain the differences in ε L . The assumption of plants being well adapted to nutrient and water stress might not be true for the restored Hor site, where parts of the vegetation still consist of species which are not typical for wetlands. This might explain the low productivity at that site, but could only be covered by a model if site specific plant responses to high water levels would be applied. Additionally, ε L is known to be species specific (Sinclair and Horie, 1989;Reich et al., 1998;Wohlfahrt et al., 1999).
Radiation-use efficiency is an important parameter needing site specific calibration. If common values were used for ε L , p ck and m retain , mean GPP would be underestimated by a factor of 2.4 (FsB) or overestimated by a factor of 3 (Lom). If site specific values were used for p ck and m retain , the discrepancy would be even higher. However ε L can easily be fitted if either GPP, biomass or LAI is known.

The control of decomposition and plant respiration by soil temperature
The annual R eco , which was dominated by summer R eco , could be described by a single temperature response function at all sites. However, it was not possible to find an equally good fit to both summer and winter R eco using the same t Q10 value. Higher t Q10 values would decrease overestimation of winter R eco especially at the southern sites FsA and FsB, but also reduce model performance for the annual R eco . Different temperature responses for different sites (e.g., Jacobs et al., 2007), seasons (e.g., Lipson et al., 2002) and temperature ranges (e.g., Lloyd and Taylor, 1994;Paul, 2001;Atkin et al., 2003) are reported in the literature. This is partly explained by multiplicative effects of several temperature sensitive processes (Davidson et al., 2006;Kirschbaum, 2006) but, still, a constant t Q10 might be a wrong assumption (Atkin et al., 2005). More sophisticated temperature responses like the Ratkowsky function (Ratkowsky et al., 1982) might improve the performance for individual sites. This might also be true for separate temperature response functions for plants and soil, as summer R eco includes autotrophic and heterotrophic respiration, while winter R eco is strongly dominated by heterotrophic respiration.

The control of decomposition by soil moisture
The activity under saturated conditions in respect to unsaturated conditions is described by p θ Satact and was strongly negatively correlated with decomposition rate k l . Patterns for p θ Satact differed between sites and variables. At all sites a minimum value of around 5 % led to the acceptable performance in annual R eco , while also quite high values did not reduce the performance except at FsB. At Lom only winter R eco was considered, as conditions were always saturated during summer. For acceptable winter R eco , p θ Satact needed to be very low. This was not true for Lom, where water in the upper soil layer partly froze in the model and led to high winter respiration.
As the soil at FsA and FsB was saturated during winter, a common lower value for p θSatact would decrease overestimation of winter fluxes. However, it would also reduce model performance at all sites and increase the site specific differences in k l (Fig. 7).
Permanently saturated soils contain less O 2 than temporally saturated ones (e.g., Kettunen et al., 1999), which affects decomposition (e.g., Reddy and Patrick, 1975;De-Laune et al., 1981;Holden et al., 2004). Therefore a lower p θ Satact would be justified for wetter sites. If k l was constant between sites and instead p θ Satact was fitted, this would lead the value of p θSatact to decrease in the order FsB > FsA > Lom > Hor > Amo (Fig. 5), which cannot be justified by the differences in water levels which increase in the order FsA < FsB < < Amo < Hor < < Lom. Therefore, a different p θ Satact cannot explain differences in soil respiration between sites. However, the amount of aerenchymatous plants, leading to soil aeration (e.g., Armstrong, 1980;Bendix et al., 1994;Grosse et al., 1996) was not taken into account. It reaches its highest coverage at FsB (90 %), followed by FsA (62 %), Hor (50 %), Lom (around 10 %) and Amo (around 6 %). Modeling water response depending on soil O 2 and redox potential, including O 2 conductance from plants, might help us to analyze the differences in decomposition rate and reduce winter overestimation. For example, in the Wetland-DNDC model, the water response function depends on redox potential: decomposition under saturated condition is reduced by a factor of only 0.6 if redox potential is high, but by a factor of 0.2 if redox potential is low (Zhang et al., 2002).

The control of decomposition by substrate
The largest differences of parameters between sites appeared for the maximum decay rate of the fast C pools, k l . Setting it to a common value would lead to an underestimation of mean R eco by a factor of 2.8 at FsB or an overestimation by a factor of 4 at Amo. Despite different temperature and water response curves being tested, k l values at FsA and FsB are substantially higher than at Amo (Figs. 5, 7). Higher t Q10 values lead to two groups of k l values: similar high ones for Lom, FsA and FsB and substantially lower ones for Hor and Amo (Fig. 7).
The partitioning into SOC pools strongly effects the differences, as can be shown by calculating decomposition rates for the total SOC (k tot ) based on k l , k h and SOC in the pools of the upper 30 cm as used in the C1 scenario (Fig. 8). However, FsB and FsA still have much higher rates than Amo. The resulting values and ranges of k tot (0.02-0.16 a −1 ) are comparable with the reported values from laboratory incubation studies of peat cores (0.03-1.66 a −1 , Moore and Dalva, 1997; 0.01-0.35 a −1 , Glatzel et al., 2004; 0.008 a −1 , Kechavarzi et al., 2010; a SOC content of 30 % was assumed for conversion from dry mass).
Lower decomposability is often associated with higher C : N ratios (e.g., Zeitz and Velty, 2002;Limpens and Berendse, 2003;Bragazza et al., 2006;Zhang et al., 2008), which might be important especially for the moss-rich Amo and Lom. Assuming a C : N ratio of 60 for the fast pools (Fig. 7, C6) leads to a decomposition rate at Lom which is close to those at FsA and FsB, while those of Hor and Amo remain substantially lower.
Low pH might be one reason for the low k l at Amo (e.g., DeLaune et al., 1981;Bergman et al., 1999). Despite being nutrient-rich and having a high pH and high biomass production, leading to large amounts of labile carbon added to the soil, k l values at Hor were very low. This might be connected to land use history and the origin of the peat from partly clayey lake sediment. Most of the labile C in the parent peat in the upper, formerly drained soil layers might have been decomposed before and therefore became stabilized.
In the current setup the slow pools were almost inert. A higher decay rate for the slow pools would result in a lower k l for sites with high C stock in the slow pools (cf. Table 3). This would decrease the differences between FsA and FsB compared to Lom and Amo, but increase the differences between FsA compared to FsB and Hor.
The highest decomposition rates occurred at sites with highest biomass production. A correlation of productivity with soil respiration was found in several comparison studies (e.g., Janssens et al., 2001;Reichstein et al., 2003). Fresh material provided by the plants might lead to higher microbial activity and priming effect (e.g., Kuzyakov, 2002;Fontaine et al., 2007). A higher plant to soil respiration ratio reduced the differences in k l between the sites and lowered winter R eco , especially at the highly productive FsA and FsB sites, but also reduced the model performance at all sites except Amo.
Vegetation at Amo and Lom consist largely of mosses which are more resistant to decomposition than vascular plants (Rudolph and Samland, 1985;Verhoeven and Toth, 1995;Limpens and Berendse, 2003;Moore et al., 2007) and might further explain the low k l value at Amo. Despite the lower biomass production, higher moss cover and higher C : N ratio compared to Hor, FsA or FsB, Lom has a relatively high decomposition rate. This can be explained by the very low dry bulk density, resulting in a low amount of C in the upper soil layers (Table 3), which are most exposed to decomposition (e.g., Fang and Moncrieff, 2005). Also, a low dry bulk density accompanies the low degree of degradation and therefore high amounts of labile carbon (e.g., Grosse-Brauckmann, 1990).
Despite the large differences in accumulated NEE (Fig. 1) between FsA and FsB, these sites have nearly identical decomposition rates. This confirms the expectations that the differences in NEE between FsA and FsB can be fully explained by the differences in the water table and the biomass and carbon stocks. Differences between sites with respect to CO 2 fluxes could be explained if, in addition to air temperature, water table and soil C and N stocks, site specific plant productivity and decomposition rates were also taken into account. Differences in nutrient availability and soil wetness could not explain the differences in plant productivity between the sites. Substrate quality and litter input, as well as pH values, are likely explanations for the differences in decomposition rates. A site specific interpretation was not needed for processes related to plant phenology, their response to temperature, allocation of new assimilates and plant respiration and litter fall rates.
The model parameters which strongly affected model performance were successfully constrained by the available long-term measurement data on NEE, partitioned into GPP and R eco , LAI and biomass, including rooting depth and root biomass at one site, water table, soil temperature and soil C and N stocks, as well as meteorological data and snow data at one site. It would have been useful if additional information about root biomass, litter fall and soil water content to validate the model performance in the corresponding processes was available for all sites. A second measurement of C and N stocks, several years after the first, as well as information about the degree of decomposition on all sites would have been very helpful to constrain decomposition rates and partitioning between SOM pools.
Some improvements in the model and its configuration were identified to obtain a better performance for simulations of GHG fluxes from treeless peatlands. Examples include separate temperature responses for plant and soil heterotrophic respiration. The static response to water saturated conditions needs to be replaced by a function that considers the change of O 2 in the soil.
The Supplement related to this article is available online at doi:10.5194/bg-12-125-2015-supplement.