McGill wetland model: evaluation of a peatland carbon simulator developed for global assessments

. We developed the McGill Wetland Model (MWM) based on the general structure of the Peatland Carbon Simulator (PCARS) and the Canadian Terrestrial Ecosystem Model. Three major changes were made to PCARS: (1) the light use efﬁciency model of photosynthesis was replaced with a biogeochemical description of photosynthesis; (2) the description of autotrophic respiration was changed to be consistent with the formulation of photosynthesis; and (3) the cohort, multilayer simple one box peat model effective one-year residence time litter its requires photon ﬂux (rain and relative and


Introduction
Over the last decade, the carbon (C) cycle in terrestrial and ocean ecosystems has been incorporated into a number of global climate simulations showing general agreement of a positive carbon cycle-climate feedback between the terrestrial biosphere and oceans and the atmosphere, but with large variations in the magnitude of the resulting CO 2 increase (Friedlingstein et al., 2006). It has been generally acknowledged that while most of the terrestrial models capture the essence of the C cycle they lack many processes and components that may be critical to a more realistic assessment (Denman et al., 2007). A recent example of a factor not included in the early coupled terrestrial C climate models that has a very high leverage on size of the positive feedback is the nitrogen cycle (e.g. Thornton et al., 2007;Zaehle et al., 2010). Additionally, land-use and land cover change, permafrost dynamics, and some critical but presently excluded ecosystems such as wetlands are believed to be important. Northern peatlands, the dominant form of wetland above ∼45 • N though they also occur in tropical regions, have not been included. This is in part because they represent < 4% of the global land F. St-Hilaire et al.: McGill wetland model: evaluation of a peatland carbon simulator surface (Gorham, 1995). There have been recent developments in adapting existing global ecosystem models, such as LPJ, for the inclusion of northern wetlands and permafrost ). These efforts hold promise but incorporating wetland hydrology and particularly organic soils remains a significant challenge (Frolking et al., 2009).
While the present day net primary production (NPP) of northern peatlands may represent < 1% of total terrestrial NPP, the amount of organic C stored in peatlands is very large relative to any other terrestrial biome or ecosystem -i.e. between ∼250 and 600 Pg C, or 10 to 20% (Gorham, 1991;Turunen et al., 2002;Tarnocai et al., 2009) of a ∼2300 Pg C total (Denman et al., 2007). The maintenance of this store of C is in large part a function of the moisture conditions of peatlands. If moisture were to change due to climate change, it is expected that the C uptake or release and methane (CH 4 ) emissions increase or decrease resulting from wetter or dryer conditions respectively (Moore et al., 1998). A change in stored C by 5% could represent 12 to 25 Pg C. Unfortunately, unlike forested and grassland ecosystem biogeochemistry models, there has been little effort in developing models of peatland biogeochemistry that are suitable for use in climate simulations (e.g. Frolking et al., 2002). In this paper we develop a model, the McGill Wetland Model (MWM), based on the general Peatland Carbon Model (PCARS; Frolking et al., 2002), but that has the same general structural and functional components as the Canadian Terrestrial Ecosystem Model (CTEM; Arora, 2003;Arora andBoer, 2005a,b, 2006), the terrestrial C model developed for inclusion in the Canadian Centre for Climate and Model Analysis (CCCma) coupled general circulation model. Eventually a MWM-like model would be incorporated into CTEM if the general climate models are sufficient to support the hydrological needs of wetland simulation in climate change scenarios.
Peat is the remains of partially decomposed plants and it accumulates because the NPP of a peatland exceeds decomposition, on average. Decomposition in peatlands is slow because of the persistence of anoxic conditions throughout most of the peat profile due to the saturated conditions inhibiting the diffusion of oxygen; therefore the hydrology of the ecosystem is critical to the cycling of C. In addition, many peatland plants, particularly the Sphagnum mosses that grow on the ombrotrophic (i.e. rain-fed, and/or nutrient poor peatlands) are much more resistant to decomposition than the foliar tissues of vascular plants (Moore and Basiliko, 2006). As litter is added to the peat profile the peatland surface continues to grow in height. As the litter decomposes it loses its original structure leading to a dramatic change in the pore size distribution at the long-term position of the water table. This effectively creates two layers of peat: a deep and thick anoxic zone called the catotelm and a shallow, thin, surface oxic zone called the acrotelm (Ingram, 1978). To simulate decomposition in peatlands it is essential that there be an adequate description of the hydrology of these layers of a peatland, particularly the day-to-day and seasonal variability in the position of the water table. In other work we have modified the Canadian Land Surface Scheme (CLASS) for the inclusion of organic soils and the estimation of the water table for both fen and bog type peatlands -the two dominant forms of northern peatlands Comer et al., 2000). Once the water table is known a model needs to be able to capture the differences in the rates of decomposition caused by the differences in anaerobic conditions down through the peat profile and the progressively more recalcitrant residual material that dominates at depth.
In addition to the reduction in decomposition in peatlands, a model of peatland C dynamics needs to account for the uniqueness in the plants that inhabit peatlands. Peatland vegetation is characterized by sedges, herbs, deciduous and evergreen shrubs, the latter often represented by ericaceous shrubs, mosses that are usually Sphagnum in the more nutrient poor acidic peatlands, and conifer trees, if trees are present. Most terrestrial ecosystem models can adequately represent photosynthesis and respiration for sedges using the function for grasses, and deciduous shrubs and conifers, but they lack the attributes of plant functional types that capture the behaviour of ericaceous shrubs and mosses. Mosses present a further problem, as they have no roots or vascular system.
The Canadian Terrestrial Ecosystem Model (CTEM) is representative of the general structure and function of terrestrial ecosystem models used in global couple climate simulations (Aurora, 2003). CTEM has three live C components: leaves, stem and roots; and two dead C components: litter and soil. Photosynthesis is based on the biogeochemical approach (Farquhar et al., 1980;Collatz et al., 1991Collatz et al., , 1992 with coupled photosynthesis-stomatal conductance and a description of moisture stress. Autotrophic respiration is the sum of maintenance respiration for the three live components and growth respiration. Heterotrophic respiration is the sum of respiration from a litter pool and a single soil pool, with base respiration rates modified by soil or litter temperature and moisture. To adapt PCARS closer to the structure and approach of CTEM we have: (1) replaced the light use efficiency approach for photosynthesis in PCARS with the biogeochemical approach used in CTEM and then developed the parameters for the biogeochemical model for typical peatland plants: sedges, ericaceous shrubs, mosses; (2) modified the description of autotrophic respiration to be consistent with the new formulation for photosynthesis; and (3) converted the cohort, multi-layer soil respiration model used in PCARS (the Peat Decomposition Model; Frolking et al., 2001) to a two-compartment litter and soil respiration model, where the soil (peat) is partitioned into an oxic and anoxic zone using an effective water table.
In this paper we first describe the model developments and then evaluate the performance of the MWM with the plant functional types for an ombrotrophic bog -the Mer Bleue peatland of the Fluxnet Canada and Canadian Carbon Project research networks (Lafleur et al., 2001(Lafleur et al., , 2003Roulet et al., 2007). We then examine the sensitivity of the model to changes in "key" environmental variables such as temperature and water table. We conclude with a brief discussion of how the model could be extended to other peatland types and how the MWM might be adapted for use in regional or global analyses.

McGill Wetland Model (MWM)
The MWM is composed of four C pools: two living matter pools -one moss and the other vascular plants (composed of leaves, sapwood and roots), as well as two dead matter pools -litter and peat. C enters the system through photosynthesis of vascular plants and mosses and leaves via either autotrophic respiration or heterotrophic respiration. The C allocation in roots and leaves and the simple growing degree-days approach for the seasonal phenology of vascular plants follow PCARS: a fixed maximum and minimum threshold, B max foliar and B min foliar , respectively, bound the foliar biomass of a given vascular plant and B root determines in turn the root biomass. Sapwood volume (B stem ) is a fixed parameter throughout the simulations. Moss capitulum biomass (B moss ) is also fixed and moss photosynthesises whenever environmental conditions permit. Once the vascular plant tissue and moss die they become litter and are decomposed for one year in a litter pool and then transferred to the peat C pool. At present the MWM has four plant functional types (PFTs): mosses, sedges, shrubs, and conifer trees. The details of the processes that are substantially changed from PCARS to MWM are described below.

Photosynthesis
MWM computes the photosynthesis for each PFT at an hourly time step based on the Farquhar biochemical approach (Farquhar et al., 1980;Collatz et al., 1991Collatz et al., , 1992. The computation for the non-vascular PFT is slightly different since mosses do not possess stomata. For mosses, a semi-empirical model including the effects of water content on photosynthetic capacity (Tenhunen et al., 1976) and on total conductance to CO 2 (Williams and Flanagan, 1998) replaces the stomatal conductance of vascular PFTs.
For all PFTs, net photosynthesis (A n ) is expressed as: where V c is the rate of carboxylation of Rubisco, C i is the intercellular CO 2 partial pressure, * is the CO 2 compensation point in the absence of mitochondrial respiration which is related to τ , the Rubisco enzyme specificity factor and oxygen concentration, [O 2 ], through * = 0.5 [O 2 ]/τ . R d is the dark respiration and V c is determined by the minimum of the rate of carboxylation when limited by Rubisco activity (W c ) or RuBP regeneration via electron transport (W j ). We use the standard formula for W c (not shown), where the key parameter in this description is V cmax25 the maximum velocity of Rubisco carboxylation at 25 • C. The rate of electron transport (W j ) (not shown) is described in Farquhar and von Caemmerer (1982). The key variable here is the potential electron transport rate J (Smith, 1937), which is a function of intercepted photon flux density (I ) and J max the maximum light-saturated rate of electron transport whose temperature dependency is outlined by Farquhar et al. (1980) and Lloyd et al. (1995). J max at 25 • C (J max25 ) is determined from a J max : V cmax ratio (Medlyn et al., 2002).

Conductance of vascular plant types
The canopy conductance (g c ) and boundary layer conductance (g b ) are required to obtain the C i of vascular PFTs: where C s is the canopy surface CO 2 partial pressure, C a the atmospheric CO 2 partial pressure, p is the atmospheric pressure, and the constants 1.4 and 1.6 consider the reduced diffusivity of CO 2 compared to water through the leaf surface and the canopy, respectively. C i is evaluated through iteration. A land surface scheme would provide the value of g b in a coupled regional or global simulation; in the stand-alone version g b is calculated with the Ball-Berry approach (Ball et al., 1987). The Jarvis approach (Jarvis, 1976) parameterized for peatlands is used to evaluate the canopy resistance (r c ), which is inversely proportional to g c . Soil matric potential ( ) used in the calculations of canopy conductance was evaluated individually for the catotelm and the acrotelm using the formulations of Campbell (1974) and Clapp and Hornberger (1978) and the parameters for peat suggested by Letts et al. (2000). A normalized water-content function, G(θ ), parameterized for peatland by Letts et al. (2000) modifies g c to account for the water stress factor: where θ lim is the residual soil-water content, θ l is the volumetric soil-water content and θ p is the soil porosity. The function is calculated independently for fibric and hemic peat and is weighted according to the root-biomass content in each of those layers. Shrub and sedge root biomass profiles from Moore et al. (2002) are used to estimate the weighting of β in our simulations. Volumetric soil-water content 3520 F. St-Hilaire et al.: McGill wetland model: evaluation of a peatland carbon simulator is evaluated at two depths (d) corresponding to the centre of fibric and hemic layers: where W is the water table depth, sat is the soil matric potential at saturation and b is the soil texture parameter of the peat layer as suggested by Letts et al. (2000).

Total conductance of mosses
For mosses, total conductance to CO 2 (g tc ) is used to find C i instead of stomatal conductance employed for vascular plants: Total conductance is determined from a least square regression described by Williams and Flanagan (1998) as: where f is the moss water content in units of g fresh moss/g dry moss (= m + 1, where m is g water/g dry moss). This relationship is only valid up to the maximum holding capacity of mosses ( maxcap ). Soil-water content and the capitulum interception of atmospheric water determine the water content of mosses. A function derived from the results of an experiment done by Hayward and Clymo (1982) with Sphagnum capillifolium determines the moss water content from capillary rise ( cr ) in g water/g dry moss: where mincap is the minimum interception capacity for mosses. The water content in the capitulum of mosses ( ca ) is added to the total moss water content ( m ): In turn, the intercepted water pool is affected by a loss rate, k d , due to evapotranspiration (Frolking et al., 1996): during a rain event; where t refers here to the hourly time steps, ρ water is an approximation of the rain water density, hppt is precipitation in mm h −1 , t d is the sum of the number of one-hour time steps with no precipitation. This sum is reset to zero as soon as a precipitation event occurs. If MWM were coupled to a surface climate model, Eqs. (10) and (11) would not be necessary since they would be derived directly from the latent heat flux.

Autotrophic respiration
The temperature dependency of the autotrophic respiration (AR) of mosses follows a Q 10 type relationship and is further modified by the function f to account for the moss water content effect on respiration (Fig. 2e-f;Frolking et al., 1996). A Q 10 of 2.0 Arora, 2003) along with the base rate respiration at 25 • C, R d25 , are used to calculate total dark respiration at temperature T (in • C): The autotrophic respiration of other PFTs also follows a Q 10 relationship for temperature sensitivity and is a combination of maintenance respiration of the leaves, stems, roots, and growth respiration similarly to CTEM (Arora, 2003). It is closely linked to the allocation of C in the plant.

Decomposition
Heterotrophic respiration (HR) in the C stored in peat is partitioned between oxic and anoxic respiration according to the position of an effective water table. It is assumed that the mass of peat above the effective water table decomposes under oxic rates through aerobic pathways, while peat below the water table decomposes at anoxic rates through anaerobic pathways. The effective water table depth, W eff , represents the position of the water table that is derived from the actual water table depth by adding the water distributed in the oxic layer expressed as depth and subtracting the air volume trapped in the anoxic layer. An hourly moisture profile is used to estimate the amount of water in the oxic compartment. Each compartment is characterized by either oxic or anoxic conditions with corresponding rates of respirations equal to: where k eff,o and k eff,a are termed the effective hourly mass loss rates in oxic and anoxic conditions, respectively, C o and C an are the carbon contents in the oxic and anoxic compartment, respectively. The temperature dependency of decomposition, f t , is similar to that used in PCARS  with the addition of a minimum temperature for decomposition (Clein and Schimel, 1995). We use the peat bulk density profile based on Fig. 1c in Frolking et al. (2001) to find the carbon content, which is also fractioned in the oxic and anoxic compartments accordingly with the effective water table depth: where PD is the total peat depth and frac is the biomass to carbon ratio. Peat depth requires initialization (PD 0 ) and is site specific. Fresh litter is decomposed in a separate compartment for a year using Eq. (13), with k eff replaced with an initial decomposition rate (k 0 ) for moss and for all other litter and Co is replaced with the mass of moss and all vascular plant litter, respectively. Total C content, or equivalent peat depth, is obtained by adding Eqs. (13) and (14), by subtracting from it the loss in C due to decomposition and adding to it the remaining litter from the plants after its initial year of decomposition, and finally by solving the quadratic equation for PD. Fresh litter C content is therefore not included to the peat C pool in its first year. The Peat Decomposition Model (PDM) developed by Frolking et al. (2001) is used to obtain a "representative" vertical profile of mass loss rates for bogs and fens. The profiles are built using the long-term fixed water table depths of Frolking et al. (2001) for a representative bog and fen, but the effect of anaerobic conditions on decomposition is kept as in PDM: a modifier equal to 0.1 for fens or 0.025 for bogs is used for anoxic conditions. During the initialization of the peat profile the peat temperature profile is also assumed constant. For MWM k eff,o and k eff,a are then obtained by integrating the area under the exponential mass loss curves of the profile in the oxic and anoxic layer, respectively (e.g. see Fig. 2; Frolking et al., 2001).

Site and data sets
The fluxes of CO 2 in the MWM, such as photosynthesis and respiration, are functions of environmental drivers. These drivers can either be input to the model from measurements from a specific site or can be obtained from a land surface model or global climate model, if MWM is being run in a coupled mode. The model requires hourly weather data: air and soil temperatures, water table depth, photosynthetic photon flux density, precipitation (rain and snow), wind speed, atmospheric pressure, atmospheric CO 2 concentration, relative humidity and net radiation.
For the purposes of the present study we run the MWM using 8 years of environmental measurements (1 January 1999 to 31 December 2006) from the Mer Bleue peatland, a 28 km 2 raised ombrotrophic bog near Ottawa, Canada (45 • 25 N, 75 • 40 W). We use the calendar year for our simulations. The climate of the region where Mer Bleue is located is cool-temperate with a mean annual temperature of 6.0 • C and a mean annual precipitation of 944mm for the period 1970-2000 (www.climate.weatheroffice.ec.gc.ca/ climate normals/index e.html). Hourly weather data is taken from the MB flux tower data set (http://fluxnet.ccrp.ec.gc.ca/ e about.htm).
A complementary data set containing model parameters based on studies reported in the literature serves for all sites within a range of general northern peatlands types (Table 1). There are three sets of parameters and initialization variables in Table 1. The first set is generic, meaning the parameters apply to any peatland that has the PFTs that are referred to. These are referred to under Values at 25 • C and Others (Table 1). The second set of parameters or initial values are specific to Mer Bleue. These are Site Specific initial values for vegetation biomass and peat depth. These values will be unique to each peatland that the MWM would be evaluated for but would be replaced by generic "representative" peatland values for bogs of a general ecoclimatic region for climate -ecosystem simulations. The third set of parameters, which contains the only parameter not obtained from the literature, was the V cmax25 for evergreen shrubs. We could not find any published values for peatland evergreen shrubs so we conducted a set of measurements over one summer to determine CO 2 response curves for Chamaedaphne calyculata and Ledum groenlandicum using a LICOR 6400 portable photosynthesis system in saturated light (∼1500 µmol m −2 s −1 ) and leaf chamber temperature of 25 • C. The value reported in Table 1 is the median value of over 50 measurements but we have yet to publish this work.

Results and discussion
We first assess how well MWM performed in capturing the annual and seasonal patterns and magnitude of C exchanges using the 8 years of continuous measurements from the Mer Bleue peatland. We examine the patterns of gross primary production (GPP), ecosystem respiration (ER), and net ecosystem exchange (NEE) and then examine the sensitivity of the MWM output to changes in the key environmental variables of moisture and temperature. Details on the measurement of NEE and how GPP and ER were derived from the NEE observations as well as the errors and uncertainties in the observations can be found in Lafleur et al. (2001Lafleur et al. ( , 2003 and Roulet et al. (2007). In the analysis presented below it should be noted that the uncertainty can be fairly large on GPP and ER derived from gap-filled NEE records for short time scales (hourly, daily) but the uncertainty gets much smaller for long time scales (annual) (Hagen et al., 2006). In the footprint of the Mer Bleue eddy covariance measurements mosses and shrubs were the dominant PFTs. While the MWM contains the parameters for peatland trees and sedge PFTs they were not activated in the present evaluation. Tests for these other PFTs will be done as NEE records for a wider and more varied range on peatlands types become available.

Annual patterns of simulated and measured exchange fluxes
We summed the daily gap-filled NEE from Mer Bleue to generate an annual net ecosystem productivity (NEP). Here we use the terminology for NEP as proposed by (Chapin et al., 2006): NEP is the difference between GPP and ER and equals -NEE. From the output of MWM we estimated net primary production (NPP) of the mosses and shrubs as the difference between their GPP and AR respectively. We can compare this simulated NPP with the annual estimates of NPP for Mer Bleue of Moore et al. (2002) and the range of NPP found in the literature for open bogs. Finally, MWM produces an output of total HR based on the sum of oxic decomposition of the first year litter and the peat located above the effective water table and anoxic decomposition from below the effective water table. At present we cannot do a complete analysis of net ecosystem C balance, NECB (Chapin et al., 2006), because we have not yet incorporated modules that partition the decomposition products into CO 2 and CH 4 fluxes, and net DOC export: currently, ER all goes to CO 2 . This means MWM annual ER should exceed, on average, the eddy covariance measurements of ER by ∼15 g C m −2 yr −1 based on the six year estimates of NECB . In general, the MWM simulates the magnitudes and interannual trend in annual NEP well ( Table 2). The maximum NEP underestimate was 64 g C m −2 yr −1 in 1999 and the maximum overestimate was 46 g C m −2 yr −1 in 2000. The average absolute difference between simulated and measured NEP is 39 g C m −2 yr −1 . NEP is underestimated for two of the eight years (1999,2006) and overestimated in the other years. GPP underestimation and overestimation followed the same pattern as NEP. The mean difference between observed and simulated NEP after 8 years of simulation is only 1 g C m −2 yr −1 , or < 20%. We did a Spearman's correlation analysis between simulated and observed annual GEP, ER and NEP. The results, all significant with a p >0.10, are 0.67, Table 2. Observed (Obs.), simulated (Sim.), and the difference (D) between observed and simulated ( ) annual NEP, GPP and ER for 8 years for the Mer Bleue peatland. The exchanges are expressed in g C m −2 yr −1 and the superscripts on moss and shrub GPP, ER and NPP and the HR components of ER indicate their fractional contribution. SD refers to the standard deviation. There are no direct measurements to evaluate how well MWM does in estimating the fractional components that make up total GPP and ER, but the proportions approximate what is generally expected ( Table 2). The fraction of moss GPP ranges between 0.33 and 0.39 (mean 0.36 ± 0.02) and shrub GPP from 0.61 to 0.67 (mean 0.64 ± 0.02) of the total. AR represents over 90% of ER, with shrub respiration and moss respiration comprising on average 64 ± 1% and 27 ± 2% respectively. Oxic zone decomposition contributes 96% of HR, consistent with the relative proportions of oxic and anoxic sources of CO 2 and CH 4 in the peat column experiments of Blodau et al. (2006). NPP, which is the difference between GPP and AR, displays a different pattern than the gross fluxes (Table 2). In the MWM simulation moss NPP represents a mean of 62% of total NPP (minimum and maximum of 49% and 89%), while shrubs NPP averages 38% (minimum and maximum of 11% and 51%). So while the contribution of moss and shrub to GPP and ER varies only slightly over the eight years (standard deviation [SD] of 0.02 and 0.01 g C m −2 yr −1 ) the NPP shows a much greater interannual variability (SD of 0.16 g C m −2 yr −1 ). This is due to the way MWM handles growth and maintenance respiration. In the case of moss, each year the GPP goes entirely to growing new moss, which is then assumed to die at the end of the growing season; whereas shrub has a biomass that requires significant maintenance respiration and hence a smaller fraction of GPP being translated into new biomass. MWM produces lower values of shrub NPP than expected. Measurement of the annual change in biomass in peatland shrub and moss is difficult, but the expected ranges based on a synthesis of peatland NPP studies  are 21-169 g C m −2 yr −1 for shrub above-ground NPP and 8-190 g C m −2 yr −1 for moss NPP, and 79-377 g C m −2 yr −1 for total NPP (assuming biomass is 50% C). For Mer Bleue, Moore et al. (2002) estimated above ground shrub and moss NPP in 1999 to be 80 and 85 g C m −2 yr −1 , respectively, while the MWM for the same year simulated 9 and 47 g C m −2 yr −1 , respectively. For the eight simulated years, the average of simulated above ground shrub and moss NPP were 95 g biomass m −2 yr −1 and 157 g biomass m −2 yr −1 , respectively. We believe this underestimation of shrub NPP occurs, in part, because of the range in which shrub foliar biomass is allowed to vary. We use the minimum and maximum values from PCARS , but the range could easily be greater with the water table variability observed over the 8-year evaluation period. There is, however, a dearth of empirical observations of the fractional components of total NPP in peatlands, and as far as we know, no one has reported on year-to-year variations in peatland biomass, be it aboveground or simply foliar.

Seasonal and interannual variability of simulated and measured exchange fluxes
Simulated GPP follows a strong annual cycle with maximum daily fluxes ranging from 5.0 g C m −2 d −1 to 6.0 g C m −2 d −1 during the peak growing season to zero during the coldest months (Fig. 1). Statistical analysis reveals a Willmott's index of agreement (Willmott, 1985;Comer et al., 2000) of 0.97 between the simulated and the tower fluxes with a systematic root mean square error (RMSE s ) and an unsystematic root mean square error (RMSE u ) of 0.07 g C m −2 d −1 and 0.27 g C m −2 d −1 , respectively. The low systematic error is somewhat misleading as the trend of measured versus modelled values is non-linear. There is a slight overestimation of simulated daily GPP for fluxes between 0 and ∼4 to 4.5 g C m −2 d −1 and an underestimation of observed larger fluxes (> 6 g C m −2 d −1 ) by 3 to 4.5 g C m −2 d −1 (Fig. 2). This weakness in capturing the full range of observed variability, especially the highest hourly fluxes, is not significant on an annual time scale. The tendency for MWM to underestimate the largest GPP is partly explained by the maximum threshold defined in the model for foliar biomass. However, the maximum foliar biomass should have a seasonal, not an hourly impact.
The average growing season water table depths and temperatures were ranked for the 8 years of simulation to see if there was any correlation with the average fluxes (Table 3a). There is a significant rank correlation of NEP and ER with water table, and ER with temperature (Table 3b). The standard deviation for the average temperatures is 0.79 • C and that for average water table depth is 0.06 m. According to the sensitivity analysis described below only the variation in temperatures significantly affects the fluxes. In general, GPP is greater in warmer years. However, there are exceptions to this trend. 2004 has the highest simulated GPP even though it corresponds to a relatively cold year and the lowest GPP is found in 1999, which has the warmest growing season.
Examining the inter-annual variability of cumulative simulated GPP (Table 2)    Bleue we made casual observations that indicated there was increased leaf litter fall of the evergreen shrubs. However, MWM does not allow the foliar biomass to go below a prescribed minimum value. The following year (2003) MWM grossly over-estimated GPP. Such a result would occur if the MWM carried over too much foliar biomass from the previous year. This would increase shrub photosynthesis by having more than expected leaf area to capture light and conversely increase moss photosynthesis due to a lack of shading by the shrubs. However, shrubs account for more than 65% of overall photosynthesis. Such findings underscore the importance of drought stress on the vascular plants, which was not something we initially considered an issue for peatland plants.
Yet, it appears that a year-to-year memory is needed to ensure a better description of the antecedent conditions for production in subsequent years. ER shows a strong annual cycle with maximum daily fluxes ranging between −4.2 g C m −2 d −1 and −5.2 g C m −2 d −1 during the growing season and fluxes of approximately −0.25 g C m −2 d −1 during the cold season. Simulated respiration has an agreement of 0.97 with the tower flux and a RMSE s and RMSE u of 0.14 g C m −2 d −1 and 0.23 g C m −2 d −1 , respectively. Simulated respiration is biased towards carbon loss compared to tower measurements, especially during the growing season (Fig. 3). There is a slight over-estimate of simulated ER for fluxes up to ∼ −4 to −4.5 g C m −2 d −1 , but for a small number of larger fluxes (i.e. < −6 g C m −2 d −1 ) MWM underestimates them by 1 to 3 g C m −2 d −1 (Fig. 2). While this underestimation of the flux cannot be directly attributed to a specific modelling approach in the MWM, it may suggest the need  As expected, warmer years tend to have larger ER fluxes. No correlation exists between the rankings of simulated ER and GPP fluxes. This reflects the fact that even though both fluxes are sensitive to the temperature in a similar manner, other environmental conditions also significantly affect the annual fluxes.
Simulated daily NEP shows a strong annual cycle with maximum daily uptakes ranging between 1.5 g C m −2 d −1 and 2.5 g C m −2 d −1 during the growing season and maximum ecosystem loss of around −0.25 g C m −2 d −1 during the cold season and approximately −1.0 g C m −2 d −1 during the growing season (Fig. 4). RMSE s was 0.12 g C m −2 d −1 , the RMSE u was 0.15 g C m −2 d −1 and the index of agreement 0.80 (Fig. 5). The NEP of 2004 and 2005 have the highest magnitudes while the lowest NEP occurs in 1999 and 2002. Larger NEP generally occurs in the warmer years. As mentioned earlier, daily NEP is not simulated but derived from the subtraction of ER from GPP. Therefore NEP has a tendency to underestimate the highest fluxes in a similar way to GPP. NEP also accumulates the errors propagated from both GPP and ER fluxes, generating a RMSE that represents a relative error twice as large as that for GPP and ER.

Sensitivity analysis
Sensitivity analyses were performed to assess the change in C fluxes with variations in the two main environmental parameters: water table depth and moisture supply through precipitation; and temperature, including air, surface and peat temperatures. This analysis serves two purposes. First, it gives an indication of what the key sensitivities are in the MWM and secondly, it provides some initial insights into the potential sensitivity of C cycling in northern peatlands to changes in climate. In the future we plan to use MWM coupled to a surface climate model to simulate the potential effects of climate change using the output of general climate simulations as input to the coupled wetland model. In this sensitivity analysis the structure of the ecosystem does not change due to competition among plant functional types even though the range of physical conditions imposed in the sensitivity analysis is, in some cases, well outside the range that would be considered climatic and hydrologic "niches" of the peatland plant functional plant types. The sensitivity analyses are done for the 8 years and averaged for that period (Table 4).
To fully cover the potential climatic changes, we imposed variations from the actual water table depth of −10 cm (wetter) to +30 cm (drier) in increments of 5 cm. A negative increment or a decrease in water table depth refers to a water table closer to the peat surface. The effects of the water table depth variations in moss C cycling occur through changes of moss water content, which is in turn used to calculate g tc and f θ . The changes in the shrub C cycling occur through variation in soil water content that affects the stomatal conductance. Our analysis shows that a modest decrease (increase) in water table depth results in slight decreases (increases) of both GPP and autotrophic respiration. The sensitivity of Table 4. The sensitivity of simulated GPP, autotrophic respiration (AR), NPP and oxic and anoxic heterotrophic respiration (HR) expressed in percent change relative to the baseline simulation (observed environmental variables: hwtd and airT are hourly mean water table depth and air temperature). A negative sign indicates a decrease relative to the baseline while a positive sign indicates an increase. autotrophic respiration for mosses is greater than that of GPP and therefore NPP increases (decreases) with a shallower (deeper) water table. The situation is reversed for shrubs. Consequently, the model favours shrub growth in a drier wetland and moss growth in a more humid one. HR is far more sensitive than the live plant derived fluxes to water table variations than plant functions are. Since the effective water table depth determines the partitioning between the much faster oxic decomposition rates and the slower anoxic decompositions rates, the total HR (oxic plus anoxic) increases when the water table moves deeper into the peat and decreases as the water table rises toward the peat surface. Even though the sensitivity of HR is much greater than other sensitivities, the magnitude of the fluxes derived from decomposition are relatively small, therefore the sensitivity of NEP to variations in HR is also small. Moss NPP is larger and it dictates the direction of change of NEP regardless of its low sensitivity to water table changes. In none of the simulated cases was the bog a net source of C to the atmosphere. Lafleur et al. (2005) explained the lack of an apparent relationship between water table and the observed changes in ecosystem respiration at Mer Bleue by the offset of both positive and negative factors on production and heterotrophic respiration with changes in water table.
For the temperature sensitivity analyses, we varied the mean from −2 • C to +5 • C in 1 • C increments. The analyses show that an increase (decrease) in temperature results in decrease (increase) in moss GPP and an increase (decrease) in moss AR. Moss AR is more sensitive to temperature change than GPP and therefore an increase (decrease) in temperature leads to a decrease (increase) in moss NPP. An increase (decrease) in temperature corresponds to a decrease (increase) in shrub NPP. The HR flux is equally responsive to temperature change: as temperature increases (decreases) the respiration increases (decreases). The changes in the fluxes with temperature are quite significant as temperature imposes an exponential impact upon C cycling. The Q 10 relationship used to determine the temperature sensitivity of AR and HR has a higher coefficient than the Arrhenius relationship describing that of GPP; therefore, the net effect is that NEP decreases (increases) as temperature increases (decreases). These analyses show that according to MWM ombrotrophic bogs could turn into net emitter of C to the atmosphere with a persistent rise in temperature of ∼5 • C.

Conclusion and prospects for MWM
MWM captures the primary C cycling processes in northern peatlands and simulates the C exchanges between peatlands and atmosphere within the acceptable errors, based on comparisons with measurements from the Mer Bleue ombrotrophic bog. Other major peatlands types include rich and poor fens, and both bogs and fens that support forest covers. MWM needs to be evaluated for these other peatland types before it can be applied for the regional to global assessment of the interactions between climate and general peatland carbon dynamics.
Our evaluation and sensitivity analysis identifies some areas of improvement of the MWM. The most critical problem we discovered lies in the way evergreen shrub foliar biomass is treated. It was not anticipated that a formulation for excess leaf loss due to drought stress would be needed. However, extended periods (e.g., > 30 days) with no precipitation during the growing seasons of both 2001 and 2002 resulted in extremely dry conditions at the surface of the peatland  and leaf drop from some shrubs towards the end of the summer of 2002. MWM limits the amount of foliar biomass within a specific range and currently has no capacity to shed an extra amount of litter due extended extremely dry periods. In other words, MWM lacks a function analogous to the drought stress function contained in many forest ecosystem models. Such a function would have resulted in a smaller amount of evergreen foliar biomass in the spring of 2003 and this would have reduced 2003 growing season production. Currently MWM has no interannual biomass memory. Unfortunately, our search of the literature reveals no studies reporting interannual variations in peatland vascular plant biomass.
We also suspect that the moisture content of the moss in the model does not become dry enough in years that experience drought. The supply of water to the moss is crudely modelled in the MWM. Once the water table drops below a certain depth -e.g. 20 to 30 cm, there is no significant capillary raise of water to the moss (Hayward and Clymo, 1982). Once this occurs the moss is kept moist only by atmospheric inputs and when there are extended periods with no rain we have observed the moss becomes very desiccated. In this case they may be out-competed by other plant functional types that can more easily extract water from depth (Gerdol et al., 2008;Breeuwer et al., 2009). The relationships among moss production, water level, temperature, and other environmental factors such as CO 2 concentration and nutrient availability and supply are complicated and have been systematically studied only in the past few years (e.g. Robroek et al., 2009). With more empirical studies in may soon be possible to develop a better parameterization for the structure and function of moss communities in MWM. We believe when MWM is coupled to the surface climate model we will be able to simulate plant and moss water losses much better.
The sensitivity analysis showed that MWM shrub autotrophic respiration was more sensitive to temperature than GPP, leading to a dramatic decrease in shrub NPP when temperature rises. With this decrease in shrub NPP it is very unlikely that the shrub community initialized in MWM would not be maintained. This large autotrophic respiration sensitivity to temperature is not realistic. Equation 12 is an exponential function and it would be more realistic if it reached an upper asymptote -i.e. described as a more Arrenhius type response where shrub NPP would approach zero in with increasing temperature stress. If no NPP for shrubs occurs for long time -e.g. multiple years, then the evergreen shrubs would not be sustained. This demonstrates that a dynamic vegetation structure is needed in MWM and this is a further development we wish to pursue. At present there is little understanding of the temperature -physiological response of peatland evergreen shrubs with changes in moisture and temperature (e.g. Weltzin et al., 2000).
MWM has no nitrogen cycle. Bogs are oligotrophic ecosystems but with increasing nitrogen deposition and/or increased mineralization there may be more available nitrogen. While nitrogen fertilization experiments on bogs show that there is little change in total NPP (Berendse et al., 2001), the addition of a nitrogen cycle in MWM would enable the examination of the affect of increased N deposition on bogs and/or climate induced changes in N mineralization on both plant production and peat and litter decomposition.
The MWM also needs further development to simulate the outputs of C as CH 4 and DOC. PCARS  has a crude formulation for the emission of CH 4 but it has not been widely tested. MWM does estimate anaerobic decomposition so the challenge is first estimating how much CH 4 is produced per mass of anaerobic decomposition and then emitting some of the produced CH 4 after oxidation along each of the transport pathways of diffusion, bubble flux and/or plant mediated transport. Roulet et al. (2007) and others studies conclude that DOC is a significant loss of carbon from peatlands. Some of the aerobic and anaerobic decomposition estimated in MWM has to support this net production of DOC and the simulation of this loss presents a number of challenges. Firstly, MWM will have to be coupled to a hydrological model that gives reasonable estimates of the loss of water through runoff, and secondly the partitioning of gross decomposition among CO 2 , CH 4 and net DOC export will have to be formulated to maintain continuity between the changes in C stores and fluxes. We are unaware of any studies that provide the process basis for the partitioning among the three C outputs for northern peatlands. There have been many studies of net DOC export, but none have related the export to gross DOC production or fraction of overall decomposition.
Water table depth is a key variable for peatland C cycling because it influences the spatial distribution of soil water content and subdivides the peat profile into oxic and anoxic compartments. In this stand-alone version of MWM, where there is no complementary calculations of water balance and energy balances, water table depth and soil climate are the direct inputs from field measurements. In order to investigate the response of northern peatlands to projected climate change, both water table depth and soil climate need to be simulated under the projected climate conditions. Therefore, our future plans are to couple the MWM to wetland-CLASS (Canadian Land Surface Scheme) to simulate the water table depth and soil climate. In addition, the empirical functions in this stand-alone version of MWM to simulate the moss water content will be replaced by more realistic evapotranspiration functions transferred from wetland-CLASS. After validating the coupled MWM-CLASS model against field measurements, MWM-CLASS will be ready to answer "whatif" questions and investigate how C cycling in northern peatlands may change due to projected climate change based on the IPCC emission scenarios. of the Mer Bleue observatory initially from a Natural Sciences and Engineering Research Council of Canada (NSERC) Strategic Grant, then from the Fluxnet Canada Research Network (CFCAS, NSERC and BIOCAP Canada), and currently by the Canadian Carbon Project (CFCAS). FSH and JW have also received graduate student support from the McGill University Global Environmental and Climate Change Research Centre (GEC3) and its predecessor the Centre for Climate and Global Change Research (C 2 GCR). SF has been supported by grant NSF ATM-0628399.
Edited by: A. Arneth