Simulating the Carbon Cycling of Northern Peatlands Using a Land Surface Scheme Coupled to a Wetland Carbon Model (CLASS3W-MWM)

Abstract

Northern peatlands store approximately one-third of the terrestrial soil carbon (C), although they cover only 3% of the global land mass. Northern peatlands can be subdivided into bogs and fens based on their hydrology and biogeochemistry. Peatland hydrology and biogeochemistry are tightly coupled to climate and, therefore, may be very sensitive to climate variability and change. To address the fate of the large peatland soil C storage under a future changed climate, a peatland C model, the McGill Wetland Model (MWM), was coupled to a land surface climate model (the wetland version of the Canadian Land Surface Scheme, CLASS3W), referred as CLASS3W-MWM. We evaluated the CLASS3W-MWM for a bog (Mer Bleue, located at 45.41°N, 75.48°W, in eastern Canada) and a poor fen (Degerö Stormyr, located at 64°11′N, 19°33′E, in northern Sweden).

CLASS3W-MWM captured the magnitude and direction of the present day C cycling very well for both bogs and fens. Moreover, the seasonal and interannual variability were reproduced reasonably well. Root mean square errors (RMSE) were <0.65 and the degree of agreements (d*) were >0.8 for the components of net ecosystem production (NEP) for both the Mer Bleue bog and the Degerö Stormyr fen. The performance of the coupled model for both bog and fen is similar to that of the stand-alone MWM driven by observed weather rather than simulated surface and soil climate. This modelling study suggests that northern peatlands are hydrologically and thermally conservative ecosystems. It was also shown that C cycling for bogs and fens was more sensitive to changes in air temperature than precipitation. Changes in temperature within the Intergovernmental Panel on Climate Change (IPCC) projected range switch the peatlands from a present-day C sink to a source, but projected changes in precipitation still maintain the peatlands as a C sink, although to a somewhat lesser degree. Increase in atmospheric CO₂ concentration enhances C sequestration for both bogs and fens. Our sensitivity analysis suggests that northern peatlands respond to changes in temperature, precipitation and doubled CO₂ concentration in a highly non-linear way. The sensitivity of C cycling in northern peatlands with respect to changes in air temperature, precipitation and the concentration of atmospheric CO₂ together is not a simple addition or subtraction of the sensitivity of the individual changes. Therefore, the sensitivity of a combination of changes in temperature, precipitation and doubled CO₂ concentration is very different from the sensitivity of peatlands to each environmental variable on their own. Our sensitivity analysis suggests that fens have a narrower tolerance to climate changes than bogs.

RÉSUMÉ  [Traduit par la rédaction] Les tourbières du Nord renferment approximativement le tiers du carbone se trouvant dans le sol terrestre, même si elles ne couvrent que 3% des terres du globe. On peut subdiviser les tourbières du Nord en tourbières hautes et en tourbières basses selon leur hydrologie et leur biogéochimie. L'hydrologie et la biogéochimie des tourbières sont intimement liées au climat et peuvent donc être très sensibles à la variabilité et au changement climatique. Pour étudier comment évoluera le stockage du carbone dans les grands terrains tourbeux sous un climat futur modifié, nous avons couplé un modèle de carbone de tourbière, le McGill Wetland Model (MWM), à un modèle climatique de surface terrestre (la version terres humides du CLASS3W canadien), c'est-à-dire le CLASS3W–MWM. Nous avons évalué le CLASS3W–MWM pour une tourbière haute (Mer Bleue, situé à 45,41°N, 75,48°O, dans l'est du Canada) et pour une tourbière basse ombrotrophe (Degerö Stormyr, situé à 64°11′N, 19°33′E, dans le nord de la Suède).

Le CLASS3W–MWM a très bien capturé la grandeur et la direction du recyclage actuel du carbone, tant pour les tourbières hautes que pour les tourbières basses. De plus, la variabilité saisonnière et interannuelle a été raisonnablement bien reproduire. Les écarts-types étaient <0,65 et les degrés de concordance (d*) étaient >0,8 pour les composantes de la production nette de l’écosystème tant pour la tourbière haute Mer Bleue que pour la tourbière basse Degerö Stormyr. La performance du modèle couplé pour la tourbière haute et la tourbière basse est semblable à celle du MWM autonome piloté par des conditions observées plutôt que par un climat simulé de la surface et du sol. Cette étude par modèle suggère que les tourbières du Nord sont des écosystèmes hydrologiquement et thermiquement conservatifs. Il a aussi été démontré que le recyclage du carbone pour les tourbières hautes et basses était plus sensible aux changements dans la température de l'air que dans les précipitations. Des changements de température de l'ordre de ceux projetés par le Groupe d'experts intergouvernemental sur l’évolution du climat (GIEC) font que les actuels puits de carbone que constituent les tourbières se transforment en sources, mais les changements projetés dans les précipitations maintiennent encore les tourbières comme des puits de carbone, quoique dans une moindre mesure. L'accroissement de la concentration du CO₂ atmosphérique améliore la séquestration du carbone à la fois pour les tourbières hautes et les tourbières basses. Notre analyse de sensibilité suggère que les tourbières du Nord réagissent aux changements dans la température et les précipitations et à une concentration doublée de CO₂ d'une façon fort peu linéaire. La sensibilité du recyclage du carbone dans les tourbières du Nord par rapport aux changements dans la température de l'air, les précipitations et la concentration du CO₂ atmosphérique ensemble n'est pas une simple addition ou soustraction de la sensibilité aux changements individuels. Par conséquent, la sensibilité à une combinaison de changements dans la température et les précipitations et à une concentration doublée de CO₂ est très différente de la sensibilité des tourbières à chaque variable environnementale prise seule. Notre analyse de sensibilité suggère que les tourbières basses ont une plus faible tolérance aux changements climatiques que les tourbières hautes.

Keywords:

• wetlands
• northern peatlands
• bogs
• fens
• C cycling
• McGill Wetland Model (MWM)
• Canadian Land Surface Scheme (CLASS)
• CLASS3W-MWM
• climate change

1 Introduction

Uncertainties in terrestrial carbon (C) cycling have been recognized as one of the key uncertainties in projecting the future climate (Huntingford et al., 2009). Terrestrial and oceanic C cycling has recently been included in the present global C-coupled climate model (Denman et al., 2007). Global C-coupled climate models all show a positive feedback (Cox et al., 2000; Friedlingstein et al., 2006; Denman et al., 2007), but the magnitude varies considerably. Uncertainties exist because several significant physical processes and terrestrial ecosystems have not been included in the current generation of C-coupled climate models: nitrogen (N) cycling, fire and peatlands (Denman et al., 2007). In our previous work we developed and evaluated a peatland carbon model, the McGill Wetland Model (MWM), for bogs (St-Hilaire et al., 2010; Wu et al., 2011) and fens (Wu, 2009; Wu et al., unpublished manuscript) based on the observed weather. In this paper, we coupled MWM to a land surface scheme and then evaluated this coupled model for bogs and fens. Using the coupled model we also investigated the sensitivity of two peatland forms, a bog and a fen, to changes in temperature, precipitation and atmospheric CO₂ concentration.

Northern peatlands store approximately 450 Gt C (1 Gt C = 10¹⁵ g C) (Gorham, 1991; Bridgham et al., 2006) because net primary production (NPP) in northern peatlands has been larger than decomposition for thousands of years (Clymo et al., 1998; Roulet et al., 2007; Nilsson et al., 2008). Covering only 3–6% of the global land mass (Matthews and Fung, 1987; Lehner and Doll, 2004), peatlands store approximately one-third of the global terrestrial soil C (Gorham, 1991). Wetlands are also the largest natural emitter of methane, and northern peatlands contribute about 10% of wetland emissions (Mikaloff-Fletcher et al., 2004).

The hydrology and biogeochemistry of northern peatlands are tightly coupled to climate (Charman, 2002; Belyea and Malmer, 2004; Bridgham et al., 2008; Ise et al., 2008). Northern peatlands are located at high latitudes where the greatest change is expected to occur under future climate scenarios (Meehl et al., 2007); hence, climate change could result in significant changes in C cycling in northern peatlands (Bridgham et al., 2008). Northern peatlands are generally subdivided into bogs and fens (Bridgham et al., 1998). Bogs receive water and nutrients solely from precipitation, whereas fens receive additional water and nutrients from surface and groundwater sources. The differences in water and nutrient inputs between bogs and fens result in differences in soil chemistry, hydrological regime, nutrient availability and vegetation structure (Bridgham et al., 1998; Chapin et al., 2004). These ecological differences lead to the differences in biogeochemical processes and greenhouse gas fluxes, and bogs and fens should respond to climate change differently (Bridgham et al., 2008). Field manipulations suggest that fens are more sensitive to climate change than bogs, especially to warming temperature and water table drawdown (Bridgham et al., 2008).

To examine the sensitivity of the C cycling in northern peatlands in response to a climate change, we developed a process-based C cycling model, the MWM in which C cycling is a function of surface and soil climate, and vegetation composition (Wu, 2009; St-Hilaire et al., 2010; Wu et al., unpublished manuscript). The MWM can simulate the vertical CO₂ exchanges between peatlands and the atmosphere. To examine the C cycling under the future climate, the soil climate, including soil temperatures, soil moisture and water table depth (WTD), need to be simulated; therefore, we coupled the MWM to a land surface climate model, the wetland version of the Canadian Land Surface Scheme (CLASS3W) (Verseghy, 1991; Verseghy et al., 1993; Comer et al., 2000; Letts et al., 2000; Ouyang et al., in press). This coupled model is referred as CLASS3W-MWM. However, we need to evaluate how well CLASS3W-MWM can simulate the present CO₂ cycling in northern peatlands.

Our objectives are (1) to evaluate the effectiveness of CLASS3W-MWM for a bog peatland and a fen peatland; and (2) to investigate the sensitivity of C cycling in northern peatlands to changes in temperature, precipitation and atmospheric CO₂ concentration. We evaluated CLASS3W-MWM at Mer Bleue (MB), an ombrotrophic bog in eastern Canada (Lafleur et al., 2001, 2003; Roulet et al., 2007), and Degerö Stormyr (DS), a minerotrophic poor fen in northern Sweden (Sagerfors, 2007; Nilsson et al., 2008; Sagerfors et al., 2008).

2 Model structure

The land surface scheme (i.e., CLASS3W) is used to simulate the surface and soil climate, including soil temperature (soilT), soil moisture and WTD, based on the surface and soil water and energy balance calculation (Fig. 1). Then, the C cycling component (i.e., MWM), based on the outputs from CLASS3W, is used to simulate the CO₂ fluxes such as gross primary production (GPP), total ecosystem respiration (TER), and net ecosystem exchange (NEE; representing net ecosystem production (NEP) in this study) between peatlands and the atmosphere. We refer to this coupled model as CLASS3W-MWM. At the present time, the coupling is unidirectional, indicating that any changes in C cycling in peatlands do not affect the simulated climate dynamics and MWM does not feed back to CLASS3W.

Fig. 1 The basic coupling structure of the coupled land surface climate model (CLASS3W) and the peatland C model (MWM) (referred as CLASS3W-MWM).

CLASS3W-MWM simulates total runoff, including surface overland runoff, subsurface and peatland lateral flow and groundwater seepage, latent, sensible and ground heat flux, temperature and water content for three soil layers, water table depth, C storage in leaves, stems and roots, one soil C pool and one fresh litter C pool that decomposes much faster than the soil C pool. The inputs required are incoming shortwave and longwave radiation, precipitation (P), air temperature (T_a ), specific humidity, wind speed, atmospheric pressure and atmospheric CO₂ concentration. The current C outputs for C fluxes are GPP and autotrophic respiration (AR) for each plant functional type and oxic and anoxic decomposition. The detailed model description can be found in Verseghy (1991), Verseghy et al. (1993), Letts et al. (2000) and Ouyang et al. (in press) for the land surface model component, and Wu (2009), St-Hilaire et al. (2010), Wu et al. (2011) and Wu et al. (unpublished manuscript) for the C cycling model component.

The soil hydrology of CLASS3W-MWM was outlined by Verseghy (1991). Parameterization of hydraulic properties for peatland soil (or organic soil) was introduced by Letts et al. (2000). A water content form of Richard's equation (Jury and Horton, 2004), is solved using classic soil-water characteristic curves, and unsaturated hydraulic conductivity is a function of soil water content (Campbell, 1974). A simple hydrological scheme was added to convert the distribution of soil moisture to a WTD based on the saturated soil volume and specific yield and retention (Letts et al., 2000). The original wetland version of CLASS works well for fens (Comer et al., 2000; Letts et al., 2000), but it does not simulate the water table depth for bogs well. Therefore, a lateral outflow was added for bogs to update the soil water dynamics in order to simulate the water table depth properly (Ouyang et al., in press).

In CLASS3W-MWM, four plant functional types (PFT), trees, shrubs, sedges, and mosses, are used to represent the vegetation communities in northern peatlands. The details of the carbon model, MWM, are provided in Wu (2009), St-Hilaire et al. (2010), Wu et al. (2011) and Wu et al. (unpublished manuscript). In summary, C enters a peatland ecosystem through photosynthesis by vascular plants, mosses and trees, if trees are present. This system loses C through AR of plants, decomposition of fresh litter, and soil decomposition (i.e., heterotrophic respiration (HR)). Photosynthesis in vascular plants is simulated for each PFT by a Farquhar biochemical approach (Farquhar et al., 1980; Farquhar and von Caemmerer, 1982), coupled with canopy-level stomata conductance (Jarvis, 1976; Ball et al., 1987) to solve an unconstrained PFT-specific GPP (Arora and Boer, 2003; Wu, 2009; St-Hilaire et al., 2010). For non-vascular plants (i.e., mosses) a semi-empirical function including the effect of the water content of moss on photosynthetic capacity (Tenhunen et al., 1976) and on total conductance of CO₂ (Williams and Flanagan, 1998) is used to replace the stomata conductance. For a fen peatland, this unconstrained GPP is constrained by a non-linear function of WTD to compute the actual GPP rate (Wu, 2009; Wu et al., unpublished manuscript). The fixed C is allocated to leaves, stems and roots, following Frolking et al. (2002). A simple degree-day approach is used for the vegetation phenology (Frolking et al., 2002; St-Hilaire et al., 2010). AR is composed of maintenance and growth respiration; HR is the sum of fresh litter decomposition and peat soil decomposition. Peat soil decomposition is calculated based on two soil carbon compartments, specifically, oxic C storage above the WTD and anoxic C storage below the WTD. For a bog peatland, the oxic and anoxic zone is subdivided by an effective WTD based on the actual WTD (St-Hilaire et al., 2010), whereas for a fen, because the roots of sedge may penetrate below the WTD, the effective WTD is calculated from a non-linear function of sedge NPP (Wu 2009; Wu et al., unpublished manuscript).

3 Evaluation study sites and datasets

a Evaluation Study Sites

The evaluation peatlands are Mer Bleue (MB), an ombrotrophic bog, and Degerö Stormyr (DS), a minerotrophic poor fen. Continuous measurements of C, water and energy exchange between the peatland surface and the atmosphere have been taken since May 1998 at MB (Lafleur et al., 2003; Roulet et al., 2007) and since 2001 at DS (Sagerfors, 2007; Sagerfors et al., 2008).

1 Mer Bleue bog

Mer Bleue is located 10 km east of Ottawa, Ontario, Canada (45.41°N, 75.48°W, 69 m above mean sea level), in a postglacial channel, and it began to accumulate organic matter 9000 years before present and changed to a bog phase between 7100 and 6800 years ago (Auer, 1930; Mott and Camfield, 1969). The climate of the region is classified as cool continental temperate, with a 30-year (1971–2000) mean annual temperature of 6.0°C and a total annual precipitation of 943 mm, 235 mm of which falls as snow between December and March (Environment Canada, 2009) (Table 1). Total precipitation is relatively evenly distributed among the 12 months, with a minimum of 60 mm in February and a maximum of 90 mm in July. The warmest month is July and the coldest is January with 30-year average monthly temperatures of 20.9°C and −10.8°C, respectively. The mean monthly temperature is >0°C from April to October and >10°C from May to September. The 1971–2000 average number of days with minimum temperature >0°C was 206 days.

Table 1. The detailed physical characteristics and mean climate conditions for the two evaluation peatlands.

Characterstics	Bog (MB)	Fen (DS)
Name	Mer Bleue	Degerö Stormyr
Location	45.41°N, 75.48°W, eastern Canada	64°11′N, 19°33′E, northern Sweden
Climate	Cool continental temperate	Cold temperate humid
Peatland classification	Ombrotrophic bog	Minerotrophic poor fen
Nutrients	Poor	Rich
Water sources	Precipitation	Precipitation, surface runoff and groundwater
Size	28 km^2	6.5 km^2
Depth of peat	5–6 m	3–4 m
Basal age	approximately 9000 years	approximately 8000 years
Dominant plants	Moss and shrubs	Moss and sedges
Surface topography	69 m above mean sea level	270 m above mean sea level
Mean annual Ta (°C)	6.0	1.2
Monthly max. Ta (°C)	20.9	14.7
Monthly min. Ta (°C)	−10.8	−12.4
Days above 0°C	206.0	150.0
Mean annual P (mm)	943.0	523.0
Precipitation as snow (mm)	235.0	233.0
Monthly max. P (mm)	90.0	65.0
Monthly min. P (mm)	60.0	25.0
Note: MB denotes Mer Bleue and DS denotes Degerö Stormyr.

Mer Bleue is a 28 km² ombrotrophic bog (Moore et al., 2002; Lafleur et al., 2003; Roulet et al., 2007) (Table 1). The measurement site is located in the northernmost arm of MB. It is a slightly domed bog with a peat depth of 5–6 m in the centre. The dominant evergreen (Chamaedaphne calyculata, Ledum groenlandicum, Kalmia angustifolia), and deciduous (Vaccinium myrtilloides) shrubs have an average height of 20–30 cm and a leaf area index (LAI) of approximately 1.3 (one-sided) (Moore et al., 2002; Bubier et al., 2006). Sedges (Eriophorum vaginatum) consitute a sparse cover and a few small trees (Picea mariana, Larix laricina, Betula populifolia) are present on hummocks. The grounds are covered by Sphagnum mosses (S. capillifolium, S. magellanicum, S. Fuscum). However, the peatland area in this study, where the eddy covariance tower and other measurement instruments for the surface and soil climate were set up, is only covered by mosses and shrubs. The total biomass above ground for vascular species, measured in 1999 and 2004, averaged 356 ± 100 g m⁻² (Moore et al., 2002) and 433 g m⁻² (Bubier et al., 2006), respectively, and the total biomass below ground, measured in 1999, was 1820 ± 660 g m⁻². Sphagnum capitulum biomass was 144 ± 30 g m⁻² in 1999 (Moore et al., 2002) and 158 g m⁻² in 2004 (Bubier et al., 2006). More details on the vegetation communities and their characteristics can be found in Moore et al. (2002) and Bubier et al. (2006).

2 Degerö Stormyr fen

Degerö Stormyr (DS), located in northern Sweden (64°11′N, 19°33′E), is an acid, oligotrophic, minerogenic, mixed mire system covering 6.5 km² (Sagerfors, 2007; Sagerfors et al., 2008) (Table 1). The mire is situated between two rivers and consists of interconnected smaller mires divided by islets and ridges of glacial till. The depth of the peat is generally between 3 and 4 m, but depths up to 8 m have been observed. The deepest peat corresponds to an age of approximately 8000 years. The climate of the region, defined as cold temperate humid, has a 30-year average annual precipitation and temperature of 523 mm and 1.2°C, respectively, and average monthly temperatures in July and January of 14.7°C and −12.4°C, respectively (Alexandersson et al., 1991; Sagerfors 2007, 2008) (Table 1). Daily mean temperature was >5°C (Moren and Perttu, 1994) for approximately 153 days for the 2001–06 study period (Sagerfors, 2007; Sagerfors et al., 2008).

The vascular plant community is dominated by Eriophorum vaginatum L., Vaccinium oxycoccos L., Andromeda polifolia L., Rubus chamaemorus L. with both Carex limosa L. and Schezeria palustris L. occurring more sparsely. Carex rostrata L. is found within the margins of this fen where there is a direct minerogenic water inflow. The surface vegetation in the wet carpets is dominated by Sphagnum majus Russ. C. Jens, and the lawns by S. balticum Russ. C. Jens. and S. lindbergii Schimp, whereas S. fuscum Schimp. Klinggr. and S. rubellum Wils. are dominant on the hummocks.

b Data Processing for the CLASS3W-MWM Drivers

The eddy-covariance of CO₂ fluxes was measured at both MB and DS along with a variety of additional environmental variables that were used as inputs for CLASS3W-MWM (Lafleur et al., 2003; Roulet et al., 2007; Sagerfors et al., 2008). These were incoming shortwave and longwave radiation, air temperature (T_a ), peatland surface temperature, precipitation, wind speed, specific humidity, atmospheric pressure and atmospheric CO₂ concentration. For CLASS3W-MWM, these inputs were measured every 30 minutes. We used the gap-filled field measurements from MB for 1999–2006 and from DS for 2001–06 as the input drivers to CLASS3W-MWM. CLASS3W-MWM was parameterized and initialized for MB and DS separately (Table 2). To represent the site-specific biogeochemical processes fully, several site-specific parameters were required: peat depth, foliar biomass for mosses, sedges and shrubs, vegetation root biomass, the maximum velocity of Rubisco carboxylation at 25°C (V _cmax25) for mosses, sedges and shrubs (Table 2). In addition to these site-specific parameters, a complementary dataset containing model parameters based on studies reported in the literature served for all sites within a range of general northern peatlands (see Table 1 in St-Hilaire et al. (2010)). Because of power outages, instrument failures and other conditions that did not satisfy the measurement requirements, gaps existed in the measurement data. Gaps for the half-hourly measurement of environmental drivers were simply filled by the average of the measurements taken over the previous 5 hours.

Table 2. Parameters and initial values used in CLASS3W-MWM for MB and DS.

Parameter	Bog (MB)	Fen (DS)	Units	Description	Reference/Source
Values for Model Parameters at 25°C
J max:V max	1.67	1.67	—	ratio	(Medlyn et al., 2002)
Mosses
V cmax25	6	6.5	μmol m^−2 s^−1	max carboxylation rate (spring)
	14	14	μmol m^−2 s^−1	max carboxylation rate (summer)
	7	7.5	μmol m^−2 s^−1	max carboxylation rate (autumn)
R d25	0.946	0.946	μmol m^−2 s^−1	dark respiration rate	(Harley et al., 1986)
Shrubs/Sedges
V cmax25	17	20	μmol m^−2 s^−1	max carboxylation rate
Initial site specific parameters
frac	48.7	48.7	%	biomass to carbon ratio
B moss	140	150	g biomass m^−2	moss capitula biomass
B min foliar	175	90	g biomass m^−2	min shrub/sedge foliar biomass
B max foliar	350	180	g biomass m^−2	max shrub/sedge foliar biomass
LAI_moss	1.35	1.05	m^2/m^2	initial moss leaf area index
LAI_shrub/sedge	0.9	0.5	m^2/m^2	initial shrub/sedge leaf area index
Root_shrub/sedge	500	450	g biomass m^−2	Initial shrub/sedge root biomass
Litter_shrub/sedge	150	250	g biomass m^−2	Initial shrub/sedge litter biomass
B stem	0.003		m^3 m^−2	shrub sapwood volume
PD0	6	4	m	initial peat depth
k o	0.05/0.2	0.05/0.2	yr^−1	initial decomposition rate for moss/shrub(sedge)	(Frolking et al., 2001)
Initial values for organic soil
Z	0.20/0.75/6.0	0.20/0.75/4.0	m	Depth to bottom of 3 soil layers	(Ouyang et al., in press)
θp	0.93/0.88/0.83	0.93/0.88/0.83		Porosity for 3 soil layers	(Letts et al., 2000)
θr	0.175/0.530/0.705	0.175/0.530/0.705		Specific retention for 3 soil layers	(Letts et al., 2000)
θm	0.04/0.15/0.22	0.04/0.15/0.22		Residual soil water content for 3 soil layers	(Letts et al., 2000)
b	2.7/6.1/9.5	2.7/6.1/9.5		Soil texture parameter	(Letts et al., 2000)
Ψs	9.6 × 10^−3/9.06 × 10^−3/8.3 × 10^−3	9.6 × 10^−3/9.06 × 10^−3/8.3 × 10^−3	m	Soil matric potential at saturation	(Ouyang et al., in press)
γ	1.7 × 10^−4/2.0 × 10^−6/1. × 10^−8	1.7 × 10^−4/2.0 × 10^−6/1.0 × 10^−8	m s^−1	Hydraulic conductivity at saturation	(Ouyang et al., in press)

c Data Processing for Tower Measurement of CO₂ Fluxes

In this study, an eddy-covariance (EC) tower measurement of NEE (net ecosystem exchange, representing net ecosystem production (NEP) in this study) and the derived fluxes, GPP and TER, were used to represent the spatial average C fluxes at the ecosystem level for MB (Lafleur et al., 2003; Roulet et al., 2007), but only NEE was used for DS (Sagerfors et al., 2008), because the derived fluxes, GPP and TER, for DS were not available. The Fluxnet-Canada measurement and processing protocols were used in our flux data processing for MB (Environment Canada, 2008). The measurement and eddy covariance (EC) flux data processing details for DS can be found in Sagerfors et al. (2008).

d Model Spin-Up

The peat decomposition model (PDM) (Frolking et al., 2001) was used to initialize a “representative” vertical profile of mass loss rates for bogs and fens for CLASS3W-MWM (St-Hilaire et al., 2010). The vertical profile was developed using a long-term constant WTD of a representative bog and fen, but a modifier representing the effect of anoxic conditions on decomposition was applied in PDM (St-Hilaire et al., 2010). In CLASS3W-MWM, a different modifier was used for a bog and fen (i.e., 0.1 and 0.025, respectively). Then for CLASS3W-MWM the effective oxic and anoxic decomposition rates for a bog and fen were derived by integrating the decomposition curve of the vertical profile in oxic and anoxic zones, respectively, in the PDM. The effective oxic and anoxic decomposition rates for a bog and fen were kept constant to compute the oxic and anoxic decomposition for MB and DS, respectively, for the entire simulation period.

4 Methodology

a Evaluation Procedures

We evaluated the performance of CLASS3W-MWM by starting with a model stability test for the MB site. We drove CLASS3W-MWM for 80 years with 1999–2006 environmental data for MB. The key environmental variables to determine and affect the soil biogeochemical processes in peatland ecosystems are WTD and soilT (Charman, 2002). Moreover, they are the two commonly measured environmental variables in peatland studies. Therefore, these two simulated key variables (i.e., WTD and soilT) were also compared with the observations to evaluate the performance of CLASS3W-MWM. We then compared the simulated CO₂ fluxes with the observed ones for MB (using GPP, TER and NEP) and for DS (using only NEP).

CLASS3W-MWM was run with a 30-minute time step for the surface and soil climate components and a 60-minute time step for the C cycling components. Then the simulated hourly CO₂ fluxes were summed to obtain the daily and annual CO₂ fluxes. CLASS3W-MWM calculated GPP and AR for mosses, sedges and shrubs separately. The PFT-specific GPPs were summed to obtain the ecosystem-level GPP; the NPP was computed by subtracting AR from GPP. The HR, composed of fresh litter decomposition, oxic decomposition (oxic HR) and anoxic decomposition (anoxic HR) was subtracted from NPP to obtain NEE. All of the respiration components, including maintenance respiration, growth respiration, fresh litter decomposition, oxic and anoxic HR, were added to give TER (i.e., TER = AR + HR). In this study, we used NEE to represent NEP. We followed the convention and terminology of Chapin III et al. (2006); hence, a positive NEP indicates a C uptake by the peatland and a negative NEP indicates a C source to the atmosphere.

b Statistics

The statistical indices used in this study to evaluate model performance were the root mean square error (RMSE), that combines the systematic (RMSEs) and unsystematic RMSE (RMSEu), where RMSE² = RMSEs² + RMSEu² (Wu, 2009; St-Hilaire et al., 2010), the coefficient of determination (r ²) and Willmott's index of agreement (d*) where 0 indicates no agreement and 1 indicates perfect agreement between simulated and observed values (Willmott, 1982, 1985). These statistical indices have been commonly used in meteorological and climate research to test and evaluate models (e.g., Comer et al., 2000; Letts et al., 2000; St-Hilaire et al., 2010; Wu et al., 2011). The Willmott degree of agreement presents the performance of a model very well because the calculation is based on the direct comparison of the simulated outputs with the observations (Willmott, 1982, 1985). So it not only defines the relationship between the simulated outputs and the observations, but also describes the discrepancy between them. A good model is indicated by a high value of d* (Willmott, 1982, 1985). Moreover, a good model should have a low RMSE with a correspondingly low RMSEs and RMSEu. The unsystematic errors, represented by RMSEu, are a measure of model parameterization that can usually be corrected with changes in model parameters. The systematic errors, represented by RMSEs, are a measure of a model's potential accuracy in its current formulation. Therefore, a good model should also have a higher RMSEu than RMSEs, because a higher value of RMSEs suggests that there is no clear relationship between the simulated values and the observations; therefore, only model reformulation, rather than parameter tuning-up, should be conducted to improve the model's performance (Willmott, 1982, 1985; Comer et al., 2000; Letts et al., 2000; St-Hilaire et al., 2010; Wu et al., 2011).

c Experimental Design of the Sensitivity Analysis

We conducted a sensitivity analysis to examine how the C cycling in northern peatlands might change in response to possible changes in temperature, precipitation, atmospheric CO₂ concentration and various combinations within the range of projections (IPCC, 2007). The sensitivity tests were carried out by imposing a change in T_a between −3°C and 5°C in increments of 1°C (a negative change means colder whereas a positive change means warmer), a change in P between −30% and 30% in 10% increments (a negative change means decreased precipitation whereas a positive change means increased precipitation), doubled atmospheric CO₂ concentration ([CO₂]) and various combinations of the above: changes in T_a with a 30% increase or decrease in P (T_a  + P); changes in T_a and a 30% increase or decrease in P with a doubled [CO₂] (T_a  + P + CO₂). The changes were imposed on the half-hourly T_a , P and CO₂ concentration for 1999 to 2006 at MB and for 2001 to 2006 at DS. The simulated annual mean GPP, TER and NEP were compared with the simulated baseline GPP, TER and NEP using the observation data for the study period.

5 Results

a Model Stability Test

To test the model's stability we ran CLASS3W-MWM for 80 years using the 1999–2006 data for MB. The WTD showed interannual and seasonal variations, but ranged between 10 and 65 cm (Fig. 2a) with an 8-year repeat pattern, indicating that the model has no drift. The first layer of soil temperature (soil1_T) also showed interannual and seasonal variations, ranging between 0 and 25°C (Fig. 2b) with no drift over the 80 years. We conclude from this analysis that variations in simulated outputs result from the variability in the model inputs and not the model structure.

Fig. 2 80-year simulation of (a) water table depth (WTD) and (b) the temperature of the first soil layer (soil1_T) based on the 10 repeated simulations of 1999–2006 for Mer Bleue (MB).

We compared the simulated summer WTD and soil1_T with the observations for MB (Figs 3a and 3b) and DS (Figs 3c and 3d). The 1:1 scatterplots and their associated r ² (>0.48 for WTD and >0.73 for soil1_T) showed that CLASS3W-MWM simulated the WTD and soil1_T reasonably well, although the regression equation indicated overestimation and underestimation bias did exist. This verification established a good foundation for surface and soil climate that were used to model the CO₂ cycling for northern peatlands using CLASS3W-MWM.

Fig. 3 Comparison of CLASS3W-MWM simulated summer WTD and the temperature of the first soil layer (soil1_T) with observations for MB ((a) and (b)) and Degerö Stormyr (DS) ((c) and (d)). The regression equation and r ² are also shown. Note: The black line in the graph is the perfect 1:1 line.

b Mean Annual Estimates of CO₂ Fluxes

CLASS3W-MWM simulated the annual GPP reasonably well for MB (approximately 2% difference between simulated and observed GPP), but a year-to-year variation did exist (Table 3), ranging between −19% and 20% (the minus sign indicates underestimation whereas the positive sign indicates overestimation). The simulated annual GPP was not significantly different from the observed one (p >> 0.05). For both MB and DS, our simulation showed that shrubs (sedges) contributed a higher percentage to the total GPP than mosses (Tables 3 and 4).

Table 3. CLASS3W-MWM simulated and observed annual GPP (GPP_M, GPP_O), TER (TER_M, TER_O), NEP (NEP_M, NEP_O) for MB; and the annual NPP and the percentage contribution of each simulated component. The suffix “O” represents observed fluxes and “M” represents simulated fluxes. The 8-year average (Ave) is also shown. For comparison, the stand-alone MWM simulated fluxes (GPP_MO, TER_MO, NEP_MO), based on observed soil temperature and WTD, are also presented (reproduced from St-Hilaire et al. (2010)). The anomalies of annual average air temperature (ΔT) and precipitation (ΔP), relative to the 8-year average value (Ave), are shown.

Year	GPP					TER							NEP			NPP
	GPP_O	GPP_M	GPP_MO	Moss	Shrub	TER_O	TER_M	TER_MO	AR	MossAR	ShrubAR	HR	NEP_O	NEP_M	NEP_MO	NPP_M	Moss	Shrub	ΔT	ΔP
	g C m^−2 yr^−1	g C m^−2 yr^−1	g C m^−2 yr^−1	%	%	g C m^−2 yr^−1	g C m^−2 yr^−1	g C m^−2 yr^−1	%	%	%	%	g C m^−2 yr^−1	g C m^−2 yr^−1	g C m^−2 yr^−1	g C m^−2 yr^−1	%	%	°C	mm
1999	646.2	607.0	624.0	42.5	57.5	−581.7	−466.3	−623.0	73.8	25.1	48.7	26.2	64.6	140.7	1.0	262.9	53.7	46.3	0.6	6.9
2000	463.2	568.7	628.0	45.1	54.9	−431.3	−475.9	−550.0	71.2	28.1	43.1	28.8	31.9	92.8	78.0	229.8	53.5	46.5	−0.9	−129.9
2001	543.3	638.5	662.0	42.5	57.5	−541.3	−496.4	−623.0	71.4	24.0	47.3	28.6	2.0	142.1	39.0	284.2	53.4	46.6	0.6	−200.1
2002	511.0	595.3	647.0	43.6	56.4	−498.1	−501.4	−612.0	70.1	25.9	44.2	29.9	12.9	93.9	35.0	244.0	53.1	46.9	0.2	−54.3
2003	494.6	592.2	667.0	40.8	59.2	−480.2	−521.9	−582.0	73.3	26.9	46.3	26.7	14.4	70.4	85.0	209.9	48.2	51.8	−1.0	−11.1
2004	683.0	619.3	713.0	41.2	58.8	−567.8	−517.9	−580.0	75.5	28.3	47.1	24.5	115.2	101.4	133.0	228.5	47.5	52.5	−0.9	167.1
2005	688.3	632.3	710.0	38.7	61.3	−597.9	−566.1	−609.0	75.9	27.6	48.3	24.1	90.4	66.2	101.0	202.7	43.6	56.4	0.1	19.5
2006	772.2	628.2	704.0	39.8	60.2	−624.7	−556.8	−604.0	76.9	28.8	48.0	23.1	147.5	71.4	99.0	200.2	44.8	55.2	1.0	201.9
Ave	600.2	610.2	669.4	41.8	58.2	−540.4	−512.8	−597.9	73.5	26.9	46.6	26.5	59.9	97.4	71.4	232.8	49.7	50.3	6.2	660.9
Note: For air temperature anomalies, negative indicates colder and positive warmer temperatures. For precipitation anomalies, negative indicates drier and positive indicates wetter conditions.

The TER was simulated a little better than GPP with a range between 5 and 10% for MB (Table 3). For both MB and DS (Tables 3 and 4), of the components of the simulated TER, AR contributed >60%, of which the AR of shrubs (sedges) contributed >40%. Oxic decomposition contributed eight times more to TER than anoxic decomposition.

Table 4. CLASS3W-MWM simulated GPP (GPP_M), TER (TER_M), NEP (NEP_M), NPP (NPP_M) and observed NEP (NEP_O), and the percentage contribution of each simulated component for DS. The suffix “O” represents observed fluxes and “M” represents simulated fluxes. The 6-year average (Ave) is also shown. For comparison, the stand-alone MWM simulated fluxes (GPP_MO, TER_MO, NEP_MO), based on observed soil temperature and WTD, are also presented (Wu et al., unpublished manuscript). The annual average air temperature anomalies (ΔT) and precipitation anomalies (ΔP), relative to the 6-year average value (Ave), are shown.

Year	GPP				TER						NEP			NPP
	GPP_M	GPP_MO	Moss	Sedge	TER_M	TER_MO	AR	MossAR	SedgeAR	HR	NEP_O	NEP_M	NEP_MO	NPP_M	Moss	Sedge	ΔT	ΔP
	g C m^−2 yr^−1	g C m^−2 yr^−1	%	%	g C m^−2 yr^−1	g C m^−2 yr^−1	%	%	%	%	g C m^−2 yr^−1	g C m^−2 yr^−1	g C m^−2 yr^−1	g C m^−2 yr^−1	%	%	°C	mm
2001	443.2	394.0	36.8	63.2	357.6	344.0	73.2	18.6	54.6	26.8	47.9	85.6	50.0	181.4	53.1	46.9	0.1	162.8
2002	494.8	593.0	29.9	70.1	467.5	532.0	68.5	17.5	51.0	31.5	60.6	27.3	61.0	174.4	37.9	62.1	0.7	58.0
2003	370.6	455.0	35.1	64.9	260.6	405.0	81.9	28.9	53.0	18.1	59.1	110.0	50.0	157.1	34.8	65.2	0.0	105.2
2004	245.1	353.0	35.9	64.1	236.8	313.0	78.8	20.5	58.3	21.2	54.9	8.3	40.0	58.4	67.6	32.4	−2.6	−14.0
2005	490.5	477.0	29.9	70.1	458.0	427.0	71.2	18.0	53.2	28.8	47.5	32.5	50.0	164.2	39.1	60.9	1.3	−60.8
2006	381.9	419.0	35.3	64.7	354.7	404.0	63.7	15.9	47.8	36.3	17.3	27.2	15.0	155.9	50.4	49.6	0.7	−251.2
Ave	404.3	449.0	33.8	66.2	355.9	404.0	72.9	19.9	53.0	27.1	47.9	48.5	44.0	148.6	47.2	52.8	2.3	568.0
Note: For air temperature anomalies, negative indicates colder and positive indicates warmer temperatures. For precipitation anomalies, negative indicates drier and positive indicates wetter conditions.

The difference between the two relatively large fluxes, GPP-TER, gives the NEP. The observed NEP represented 0.009–0.20 of GPP or TER and the discrepancy between the observed and simulated NEP was within the uncertainty of both GPP and TER (Tables 3 and 4). Over the study period, the difference between observed and simulated NEP was within the standard deviation of the measurement reported in Roulet et al. (2007) for MB, Sagerfors (2007) and Sagerfors et al. (2008) for DS.

The year-to-year variation of simulated NPP (Tables 3 and 4) was within the observed range reported by Moore et al. (2002) for a bog and Nilsson and Bohlin (2001) for DS. For example, for a bog Moore et al. (2002) reported a shrub above-ground NPP of 21–169 g C m⁻² yr⁻¹; a moss NPP of 8–190 g C m⁻² yr⁻¹; and a total NPP (assuming biomass is 50% C) of 79–377 g C m⁻² yr⁻¹.

Further, the CLASS3W-MWM simulated C fluxes, including GPP, TER and NEP, for the study period (i.e., 1999–2006 for MB and 2001–06 for DS) were comparable to the ones simulated by the stand-alone MWM, where the observed soilT and WTD were used as model inputs (Tables 3 and 4) (St-Hilaire et al., 2010; Wu et al., unpublished manuscript). For example, at MB, for 1999–2006, NEP simulated by CLASS3W-MWM ranged from 66 to 141 g C m⁻² yr⁻¹ with an average of 97 g C m⁻² yr⁻¹, while NEP simulated by the stand-alone MWM was 1–133 g C m⁻² yr⁻¹ with an average of 71 g C m⁻² yr⁻¹. At DS, NEP simulated by CLASS3W-MWM ranged from 8 to 110 g C m⁻² yr⁻¹ with an average of 48 g C m⁻² yr⁻¹, while NEP simulated by the stand-alone MWM was 15–61 g C m⁻² yr⁻¹ with an average of 44 g C m⁻² yr⁻¹. Moreover, the discrepancy between the simulated and observed C fluxes did not have a specific relationship with warmer or drier climate conditions (Tables 3 and 4).

c Seasonal and Interannual Variability of CO₂ Fluxes

CLASS3W-MWM captured the interannual and seasonal variability of GPP, TER, NEP at MB and NEP at DS reasonably well (Figs 4 and 5), although it underestimated the summer peak GPP and TER at MB. CLASS3W-MWM underestimated the summer peak GPP and TER by 1–2 g C m⁻² d⁻¹ at MB. The discrepancy between the simulated and observed NEP was between −1 and +1 g C m⁻² d⁻¹, where the positive value represents overestimation and the negative value represents underestimation. In general, the maximum discrepancies occurred in the summer when both vegetation photosynthesis and soil respiration are much higher than other seasons.

Fig. 4 Temporal trend and comparison of daily observed and CLASS3W-MWM simulated GPP (top panel), TER (middle panel) and NEP (bottom panel) for 1999–2006 for MB (the filled black circles show the observed CO₂ fluxes, and the open red circles show the CLASS3W-MWM simulated CO₂ fluxes).

Fig. 5 Temporal trend and comparison of daily observed and CLASS3W-MWM simulated NEP for 2001–06 for DS (the filled black circles show the observed NEP, and the open red circles show the CLASS3W-MWM simulated NEP).

Further, we have also investigated the limits on the shrub foliar biomass that was suggested by St-Hilaire et al. (2010) as a potential source of this discrepancy. Our results showed that relaxation of the limits on the shrub foliar biomass did not improve the performance significantly. When the maximum shrub foliar biomass was increased by 100%, the coefficient of determination (r ²) for 1999–2006 between the observed and simulated GPP was 079–0.94, TER 0.62–0.92, and NEP 0.54–0.73. Comparably, with the original value of maximum shrub foliar biomass, r ² for 1999–2006 for GPP was 0.87–0.96, for TER 0.67–0.93, and for NEP 0.62–0.85. Increasing shrub biomass reduced the moss GPP because more light was intercepted by shrubs; thus, less light was available for moss photosynthesis. Therefore, GPP is the same even when the maximum shrub foliar biomass was increased by 100%.

The statistical indices (r ² > 0.40 for NEP, r ² > 0.85 for GPP and TER, and d* > 0.80 for GPP, TER and NEP) indicated that CLASS3W-MWM can simulate the CO₂ exchange of a bog and fen reasonably well (Figs 6 and 7), although the regression showed that bias did exist (Figs 6 and 7). The statistical performance of NEP was much lower than that of GPP or TER, but still well within acceptable limits. Moreover, in all cases, the systematic components of RMSE were much less than the unsystematic ones, indicating that better parameters rather than changes in model structure could improve the model's performance (Comer et al., 2000; Letts et al., 2000).

Fig. 6 The 1:1 scatterplots of observed and CLASS3W-MWM simulated daily GPP and TER for MB for 1999–2006. The solid black line is the perfect 1:1 line. The filled circles show the GPP, and the open circles show the TER. The regression equation r ², d*, and the RMSE are also shown.

Fig. 7 The 1:1 scatterplots of observed and CLASS3W-MWM simulated daily NEP for MB for 1999–2006 (top panel) and DS for 2001–06 (bottom panel). The solid black line is the perfect 1:1 line. The regression equation r ², d*, and the RMSE are also shown.

d Sensitivity Analysis

Total GPP and TER decreased as air temperature (T_a ) decreased, since each component of GPP and TER decreased (Tables 5a and 5b). Because the decrease in GPP was greater than the decrease in TER, NEP decreased slightly as T_a decreased, but remained positive (i.e., net C uptake). With an increased T_a , the total GPP for MB increased slightly in spite of a decreased moss GPP as a result of an increase in shrub GPP, whereas the total GPP for DS decreased because both moss and sedge GPP decreased. NEP decreased with an increased T_a as a result of a greater increase in TER. With an increase of 5°C for MB, and 3°C for DS, NEP became negative (i.e., a source of CO₂ to the atmosphere).

Table 5. Sensitivity of CLASS3W-MWM simulated GPP, AR, NPP , NEP, TER, OxicHR and AnoxicHR to changes in air temperature (T_a ), precipitation (P), atmospheric CO₂ concentration ([CO₂]) and combinations of these, expressed as a relative change in percentage of the base line fluxes, for Mer Bleue (MB) (a) and Degerö Stormyr (DS) (b). Negative percentage changes indicate a decrease, and positive percentage changes indicate an increase relative to the base line value. The percentage change is calculated by “(new value-base line)/base line × 100”. The baseline C fluxes were calculated based on CLASS3W-MWM using present-day climate inputs. Changes in T_a were imposed by its absolute changed values, whereas changes in precipitation were imposed by its relative changed percentage. Positive changes in T_a indicate warmer temperatures, and negative changes in T_a indicate cooler temperatures. Negative changes in precipitation indicate decreased precipitation, and positive changes indicate increased precipitation. Only doubled atmospheric CO₂ concentration was applied (2 × CO₂). Care should be taken that the changes were evenly distributed in the study period. The last three columns show the absolute C fluxes for GPP, TER and NEP when those changes were imposed.

	Base line (g C m^−2 yr^−1)	Moss (%)			Shrub (%)			HR (%)		Total (%)			Total	gC m^−2 yr^−1
		GPP	AR	NPP	GPP	AR	NPP	Oxic HR	Anoxic HR	NEE	TER	GPP	NEE	TER	GPP
		254.7	138.0	116.7	355.5	239.5	116.0	122.9	12.5	97.4	512.8	610.2	97.4	512.8	610.2
(a) For MB
	Ta (−3)	−10.8	−16.7	−3.9	−21.3	−14.2	−35.8	−19.8	−26.1	−18.9	−16.5	−16.9	78.9	428.2	507.1
	Ta (−2)	−7.1	−9.0	−4.7	−13.9	−9.1	−24.0	−12.8	−7.1	−17.3	−9.9	−11.1	80.6	462.1	542.7
	Ta (−1)	−3.0	−2.5	−3.5	−6.7	−4.5	−11.3	−6.3	−1.4	−9.5	−4.3	−5.2	88.1	490.6	578.7
T	Ta (+1)	−1.6	2.8	−6.8	6.2	6.2	6.1	5.9	7.5	−9.5	5.3	2.9	88.2	539.9	628.1
	Ta (+2)	−6.5	5.7	−21.0	12.0	13.4	9.1	9.9	12.8	−28.5	10.5	4.3	69.6	566.6	636.2
	Ta (+3)	−11.2	6.0	−31.5	16.3	23.6	1.4	16.1	22.3	−59.3	17.0	4.9	39.6	600.3	639.9
	Ta (+5)	−29.0	−2.1	−60.8	22.7	40.8	−14.7	31.1	29.3	−133.4	26.7	1.1	−32.5	649.7	617.2
	P(−30%)	−21.5	−30.6	−10.8	0.2	−1.1	2.9	11.1	12.8	−25.2	−5.8	−8.9	72.8	483.2	556.0
	P(−20%)	−7.6	−14.0	−0.050	0.1	−1.0	2.3	4.7	10.0	−4.7	−2.8	−3.1	92.8	498.3	591.1
P	P(−10%)	−1.4	−5.7	3.8	0.037	−0.4	0.9	2.4	6.4	1.6	−1.0	−0.5	99.0	507.9	606.9
	P(+10%)	0.024	4.7	−5.5	−0.025	0.3	−0.7	−1.9	−4.3	−4.4	0.8	−0.005	93.1	517.1	610.2
	P(+20%)	−0.3	7.8	−9.8	−0.041	0.8	−1.8	−3.2	−18.0	−7.5	1.3	−0.1	90.1	519.3	609.4
	P(+30%)	−0.5	9.9	−12.7	−0.1	1.2	−2.7	−4.0	−20.7	−10.8	1.8	−0.2	86.8	522.0	608.8
CO2	2 CO2	15.4	0.0	33.6	37.2	21.1	70.4	35.3	0.045	79.6	18.3	28.1	174.9	606.9	781.7
T+P	2 + (+30%)	0.5	23.6	−26.9	11.9	14.6	6.2	4.3	2.9	−30.6	14.3	7.1	67.6	586.0	653.6
	3 + (+30%)	0.1	30.0	−35.4	16.2	20.8	6.8	7.5	6.2	−44.6	19.7	9.5	53.9	614.1	668.0
	5 + (+30%)	−0.8	43.2	−52.9	22.5	39.7	−13.2	18.0	15.7	−103.9	34.9	12.7	−3.8	691.8	688.0
	2 + (−30%)	−48.0	−46.2	−50.2	12.3	16.4	3.7	66.6	18.4	−142.2	11.7	−12.9	−41.0	572.6	531.6
	3 + (−30%)	−59.5	−56.4	−63.3	17.1	23.9	3.0	79.5	20.5	−175.3	15.6	−14.9	−73.3	592.7	519.4
	5 + (−30%)	−85.5	−84.9	−86.3	24.7	47.6	−22.7	289.6	−14.6	−494.2	68.5	−21.3	−383.9	863.9	480.1
T + P +  CO2	2 + (+30%) + 2 CO2	16.8	23.6	8.9	54.8	42.5	80.2	43.2	2.9	51.2	36.6	39.0	147.3	700.7	848.0
	3 + (+30%) + 2 CO2	16.9	30.0	1.5	61.6	54.1	77.1	49.4	6.2	30.4	45.4	43.0	127.0	745.4	872.4
	5 + (+30%) + 2 CO2	16.9	43.2	−14.2	72.0	77.9	59.7	61.7	15.7	−25.8	63.2	49.0	72.3	836.9	909.2
	2 + (−30%) + 2 CO2	−39.7	−46.2	−31.9	55.0	44.4	76.9	101.0	18.4	−76.6	33.0	15.5	22.8	682.0	704.8
	3 + (−30%) + 2 CO2	−52.8	−56.4	−48.5	62.1	57.2	72.2	115.5	20.5	−120.6	39.7	14.1	−20.0	716.5	696.5
	5 + (−30%) + 2 CO2	−83.0	−84.9	−80.7	73.2	83.5	52.1	323.7	−14.6	−441.5	93.4	8.0	−332.6	991.8	659.3
	Base line (g C m^−2 yr^−1)	Moss (%)			Sedget (%)			HR (%)		Total (%)			Total	gC m^−2 yr^−1
		GPP	AR	NPP	GPP	AR	NPP	Oxic HR	Anoxic HR	NEE	TER	GPP	NEE	TER	GPP
		157.8	71.1	86.7	322.3	224.5	97.8	100.0	20.7	63.9	416.3	480.1	63.9	416.3	480.1
(b) For DS
	Ta (−3)	−7.8	−27.4	8.4	−28.6	−21.7	−44.4	−27.1	−14.3	−9.6	−23.6	−21.8	57.7	317.9	375.7
	Ta (−2)	−7.0	−18.7	2.7	−18.6	−13.9	−29.4	−15.7	−5.3	−15.3	−14.7	−14.8	54.1	355.0	409.1
	Ta (−1)	1.1	−9.6	9.9	−7.3	−4.4	−13.8	−5.4	−7.4	3.0	−5.7	−4.5	65.8	392.6	458.4
T	Ta (+1)	−1.5	12.6	−13.1	6.8	−1.5	25.7	9.9	−13.2	10.3	3.1	4.0	70.4	429.1	499.5
	Ta (+2)	−14.7	13.3	−37.6	−1.5	−2.9	1.5	32.3	−23.0	−91.9	7.3	−5.9	5.2	446.8	452.0
	Ta (+3)	−35.2	−5.2	−59.8	−17.5	−8.2	−38.7	48.3	−46.6	−200.9	4.0	−23.3	−64.4	432.7	368.3
	Ta (+5)	−48.1	−28.2	−64.5	−38.2	−3.9	−116.9	152.1	−61.2	−485.0	26.6	−41.4	−245.9	527.0	281.1
	P(−30%)	5.8	3.3	7.9	2.1	−5.4	19.4	−4.5	−16.2	52.8	−4.3	3.3	97.6	398.5	496.1
	P(−20%)	7.2	1.5	11.9	2.9	−4.0	18.7	−6.1	−16.9	59.8	−4.2	4.3	102.0	398.8	500.8
P	P(−10%)	4.9	0.7	8.3	1.4	−2.1	9.4	−3.0	−9.4	33.5	−2.2	2.5	85.2	407.0	492.3
	P(+10%)	−1.4	−0.2	−2.4	−1.6	−0.8	−3.4	−1.6	−3.6	−4.8	−1.0	−1.5	60.8	412.1	472.9
	P(+20%)	−3.6	−0.047	−6.5	−3.6	−2.4	−6.3	−3.2	−5.0	−11.8	−2.3	−3.6	56.3	406.6	462.9
	P(+30%)	−3.2	−0.027	−5.8	−2.9	−2.0	−5.1	−2.8	−3.9	−10.0	−1.9	−3.0	57.5	408.2	465.7
CO2	2 CO2	12.3	0.000	22.4	23.5	11.3	51.5	37.9	−13.8	54.4	14.5	19.8	98.6	476.7	575.3
T + P	2 + (+30%)	1.2	25.6	−18.8	14.1	6.5	31.4	21.0	−9.6	−7.3	12.5	9.8	59.2	468.1	527.3
	3 + (+30%)	−15.1	17.3	−41.6	5.5	4.7	7.3	45.1	−18.8	−109.9	15.4	−1.3	−6.3	480.2	473.9
	5 + (+30%)	−31.8	−1.0	−57.0	−6.4	5.3	−33.2	102.9	−40.0	−276.4	25.4	−14.7	−112.7	522.1	409.4
	2 + (−30%)	−40.9	−28.9	−50.8	−25.8	−24.3	−29.3	16.0	−58.1	−120.1	−17.1	−30.8	−12.8	345.2	332.4
	3 + (−30%)	−68.1	−66.6	−69.3	−44.0	−30.1	−75.9	96.8	−86.8	−333.8	−8.7	−51.9	−149.3	380.2	230.9
	5 + (−30%)	−81.8	−77.6	−85.3	−42.1	−19.3	−94.4	250.6	−78.5	−627.3	32.6	−55.2	−336.8	552.1	215.3
T + P + CO2	2 + (+30%) + 2 CO2	13.8	25.6	4.1	43.4	20.8	95.1	62.6	−24.3	61.1	29.4	33.6	102.9	538.8	641.7
	3 + (+30%) + 2 CO2	−4.1	17.3	−21.6	32.3	19.0	62.9	82.4	−29.5	−52.5	31.5	20.3	30.3	547.4	577.7
	5 + (+30%) + 2 CO2	−22.2	−1.0	−39.5	15.8	17.4	12.0	138.2	−52.6	−234.6	39.8	3.3	−86.0	581.9	495.9
	2 + (−30%) + 2 CO2	−33.5	−28.9	−37.3	−10.1	−17.0	6.0	36.6	−63.9	−78.0	−8.5	−17.8	14.1	380.8	394.9
	3 + (−30%) + 2 CO2	−64.3	−66.6	−62.3	−35.7	−25.3	−59.5	102.1	−86.7	−307.6	−4.8	−45.1	−132.6	396.2	263.6
	5 + (−30%) + 2 CO2	−80.0	−77.6	−81.9	−38.4	−16.9	−87.9	254.6	−80.2	−618.4	34.8	−52.1	−331.1	561.1	230.1
The notation is P(−30%) indicates a decrease in P of 30%; Ta (−3) indicates a decrease in Ta of 3°C; 2 CO2 inidcates doubled CO2 concentration; 2+(+30%) inidcates a combination of an increase in Ta of 2°C and an increase in P of 30%; 2+(+30%) + 2 CO2 indicates a combination of an increase in Ta of 2°C, an increase in P of 30% and doubled CO2 concentration.

Changes in precipitation (P) alone did not significantly affect the C cycling of either MB or DS, whereas peatlands were more vulnerable to decreased P than increased P (Tables 5a and 5b). Moreover, changes in P had a greater effect on mosses than shrubs or sedges. NEP increased with a slight decrease in P (e.g., a 10% decrease) but decreased in other cases. However, NEP remained positive (i.e., net C uptake) for all changes in P.

CO₂ cycling in peatlands was sensitive to changes in atmospheric CO₂ concentration ([CO₂]) (Tables 5a and 5b). Doubled [CO₂] resulted in a significant increase in each component of NEP, except moss AR and anoxic HR. Because total GPP increased more than TER, it resulted in a larger C sink and enhanced peat accumulation. However, this could result in increased dissolved organic C (DOC) export, especially for a fen (e,g., DS), which has not been considered yet in CLASS3W-MWM.

We conducted a sensitivity analysis for a combination of increased T_a (≥2°C) and increased or decreased P (±30%) (T_a  + P). Changes in T_a  + P resulted in significantly different effects in CO₂ cycling than changes in T_a or P alone (Tables 5a and 5b), although the general trends were similar to the results for changes in T_a alone. However, the effects in CO₂ cycling resulting from changes in T_a  + P were not a simple summation of the effects from changes in T_a and P individually. An increase in T_a and a 30% decrease in P resulted in both MB and DS becoming significant C sources to the atmosphere. A small increase in T_a (≤2°C for DS and ≤3°C for MB) and an increase in P still kept the peatlands functioning as a C sink, although a larger increase in T_a and an increase in P could switch the peatlands to a C source.

We also conducted a sensitivity analysis for a combination of increased T_a (≥2°C), increased or decreased P (±30%) and doubled [CO₂] (T_a  + P + CO₂) to examine the response of CO₂ cycling to a full spectrum of projected future climate scenarios (Tables 5a and 5b). The response of CO₂ cycling to T_a  + P + CO₂ was significantly different from the response to changes in T_a , P or [CO₂] alone and that of T_a  + P (Tables 5a and 5b). With a smaller increase in T_a (≤3°C for MB and ≤2°C for DS) and increased P, the effects of doubled [CO₂] offset the decreases leading to an increased C uptake. However, with a larger increase in T_a and decreased P, the C sink was reduced, and even peatlands could become a C source with a ≥3°C increase in T_a , although the magnitude was partially offset by the doubling of [CO₂].

6 Discussion

The evaluation conducted in this study by comparing the simulated C fluxes with observed ones provides us with confidence that CLASS3W-MWM can simulate the C cycling for both bogs and fens in northern peatlands reasonably well, although CLASS3W-MWM does not reproduce the summer peak GPP accurately for either bogs or fens. The statistical indices, derived from the CLASS3W-MWM simulation, are comparable with the ones from the stand-alone MWM simulation where the observed surface and soil climates were used as inputs. The values of d* for GPP, TER, and NEP at MB from CLASS3W-MWM are 0.98, 0.92, and 0.88, respectively, whereas from the stand-alone MWM are 0.97, 0.97 and 0.80, respectively (St-Hilaire et al., 2010). The values of d* for NEP at DS from CLASS3W-MWM is 0.80, whereas from the stand-alone MWM is 0.84 (Wu et al., unpublished manuscript). The RMSE for NEP at MB from CLASS3W-MWM is 0.54 g C m⁻² d⁻¹, whereas from the stand-alone MWM it is 0.23 g C m⁻² d⁻¹ (St-Hilaire et al., 2010). The RMSE for NEP at DS from CLASS3W-MWM is 0.54 g C m⁻² d⁻¹, whereas from the stand-alone MWM it is 0.49 g C m⁻² d⁻¹ (Wu et al., unpublished manuscript). Moreover, we acknowledge that the physical surface and soil climate simulated by CLASS3W-MWM are used as inputs to the simulation of C fluxes and, therefore, will translate into the errors of the C fluxes.

Our modelling exercise shows that northern peatlands are a thermally and hydrologically conservative ecosystem. Changes in air temperature do not translate directly into similar changes in peat temperature (Fig. 8). Changes in air temperature also lead to changes in soil moisture and, thus, WTD through changes in evapotranspiration. In particular, large increases in air temperature cause a dramatic drop in WTD (Fig. 8). Further, large changes in precipitation result in relatively small changes in WTD, particularly in the fen (Fig. 8).

Fig. 8 The changes in peat temperature (T: °C) and WTD (cm) relative to the baseline simulated peat temperature and WTD responding to changes in air temperature (T_a : °C) (top two panels) and changes in precipitation (P: %) (bottom two panels). Note: (1) Positive changes in peat temperature indicate an increase, and negative changes indicate a decrease; positive changes in WTD indicate an increase (or drier conditions), and negative changes indicate a decrease (or wetter conditions). (2) Different scales are used for changes in WTD for MB and DS.

Cycling of C at a bog shows greater sensitivity to increased air temperature than to decreased air temperature because increased air temperature imposes other constraints on C cycling through a possible decrease in canopy water content and soil moisture (represented by the lower WTD in Fig. 8). Moreover, respiration is an exponential function of air temperature (Frolking et al., 2002; Lafleur et al., 2005; St-Hilaire et al., 2010); therefore, changes in air temperature lead to larger changes in TER than in GPP. As a result, increased air temperature can cause a bog to become a C source to the atmosphere.

The sensitivity to air temperature of the coupled CLASS3W-MWM contrasts with the results of the stand-alone MWM for both bogs (St-Hilaire et al., 2010; Wu et al., 2011) and fens (Wu, 2009; Wu et al., unpublished manuscript). In CLASS3W-MWM, soil climate is dynamically linked with the surface climate, but changes in soil climate do not occur either consistently or simultaneously with changes in surface climate, reflected in changes in air temperature and precipitation (Fig. 8), because of the hydraulic and thermal properties of organic soil in the peatland ecosystem (Comer et al., 2000; Letts et al., 2000; Parmentier et al., 2009). For example, increases in air temperature result in increases in soil water stress (represented by increases in WTD (Fig. 8), which reduces the moss photosynthetic activity. Strack and Price (2009) and Tuittila et al. (2004) have empirically observed similar physical and biogeochemical changes with respect to changes in air temperature in the field. Our modelling exercise shows that moss GPP decreases as air temperature increases (Tables 5a and 5b). But the stand-alone MWM cannot capture this because changes in air temperature in the stand-alone MWM do not affect the soil moisture.

Moreover, our sensitivity analysis shows that the fen is more vulnerable to environmental changes than the bog. A manipulative field experiment, which studied the responses of C cycling in both bogs and fens to changes in WTD and T_a , showed a similar result (Bridgham et al., 1995, 2008). Bogs appear to tolerate a large range of environmental conditions while continuing to function reasonably normally. But fens are less resilient and tolerate a narrower range of environmental conditions and, therefore, tend to change much more dramatically than bogs with respect to the same range of environmental changes.

Our sensitivity analysis also indicates that mosses on the bog are sensitive to climate changes. Empirical studies have shown similar results through changes in moss water content (i.e., moss capitulum moisture content) (Tuittila et al., 2004; Strack and Price, 2009). It has been shown that both WTD or moss capitulum moisture content can be used to predict moss C dynamics (Strack and Price, 2009). However, it is complicated to translate changes in precipitation into changes in moss capitulum moisture content and, hence, moss photosynthetic and respiration activity. Changes in moss function have been shown to occur without apparent changes in soil water content or water table position (Robroek et al., 2009; Strack and Price, 2009). Roulet et al. (2007) speculated that the moss surface becomes effectively decoupled from the peat moisture content if the WTD is deeper than 35 cm. In CLASS3W-MWM we use both interception water from precipitation and moss capitulum water, which is a function of WTD, developed by Hayward and Clymo (1982), to represent the moss water dynamics. Because of the sensitivity of northern peatlands to moss function this is an important aspect for future study.

Our sensitivity analysis indicates that both moss and vascular plants in northern peatlands are sensitive to changes in the concentration of atmospheric CO₂. Specifically, the doubled CO₂ concentration results in increased NEP and thus enhanced peat accumulation. Field manipulations also showed similar results of elevated CO₂ enhancing the height growth of moss but not significantly influencing the growth of vascular plants (Heijmans et al., 2001). The CO₂ fertilization in northern peatlands is similar to that in other terrestrial ecosystems (Luo et al., 2008; Balshi et al., 2009); however, in CLASS3W-MWM, moss AR does not respond to changes in CO₂ concentration (Tables 5a and 5b). Moss AR is related to maximum GPP but is not dynamically linked to the actual GPP (Wu, 2009). Therefore, moss AR is simulated by a function of air temperature and moss water content only (St-Hilaire, et al., 2010; Wu, 2009). In addition, care must be taken that the growth of moss and vascular plants is not limited by nutrient availability in CLASS3W-MWM with the doubling of CO₂ concentration, although field studies have shown that nutrient availability is important to constrain the effect of elevated CO₂ (Hoosbeek et al., 2001).

The area where the northern peatlands are located will likely experience a warmer and possibly drier climate in the twenty-first century (IPCC, 2007; Meehl et al., 2007). So future climate change could have a significant impact on the hydrological and biogeochemical processes of northern peatlands (Bridgham et al., 1995, 2008; Belyea and Malmer, 2004). Our sensitivity analysis shows that changes in air temperature, precipitation and the concentration of atmospheric CO₂ together have a complex non-linear effect on the C cycling in northern peatlands. For example, decreased precipitation exacerbates the decrease in NEP resulting from increased air temperature, while increased precipitation ameliorates the decrease in NEP resulting from increased air temperature. Moreover, CO₂ fertilization can compensate for some of the decrease in NEP resulting from increased air temperature and/or decreased precipitation (Luo et al., 2008; Balshi et al., 2009). However, the sensitivity of C cycling in northern peatlands with respect to changes in air temperature, precipitation and the concentration of atmospheric CO₂ together is not a simple addition or subtraction of the sensitivity of the individual changes. Moore et al. (1998) conjectured that a non-linear relationship between hydrological and biogeochemical processes of northern peatlands introduces considerable uncertainties over exactly how northern peatlands would respond to these changes. CLASS3W-MWM could be used to carry out a factorial analysis to examine the relative role of changes in air temperature, precipitation and CO₂ concentration in maintaining northern peatlands as a C sink or changing northern peatlands into a C source to the atmosphere (e.g., Luo et al., 2008), but this is beyond the scope of the current study. Our sensitivity analysis indicates that northern peatlands are very sensitive to increased air temperature and doubling of the CO₂ concentration. These synergies should be considered when field manipulative experiments are designed.

We acknowledge that CLASS3W-MWM has a number of limitations. First, CLASS3W-MWM assumes fixed vegetation. There is no change in the structure of vegetation communities with environmental changes. However, plant communities in northern peatlands do change when they are, in particular, subjected to significant changes in moisture (Weltzin et al., 2000, 2001, 2003). Second, our sensitivity analysis shows that fens seem to be more sensitive to environmental changes than bogs, but we have assumed changes in precipitation only, and fens are also influenced by lateral exchanges in water. Therefore fens need to be put into the context of their catchment, as in the study of Sonnentag (2008) for a bog. Finally, the coupled CLASS3W-MWM has also inherited the weakness of the original MWM. One possibility for the processes not captured by the model, which may result in the simulated interannual variations in GPP, TER, and NEP not accurately following the trends of the observed ones, is that the MWM does not have an interannual biomass memory because MWM limits the amount of foliar biomass within a specific range and currently has no capacity to shed an extra amount of litter resulting from extended extremely dry periods (St-Hilaire et al., 2010; Wu et al., 2011). The other possibility is that MWM does not capture the moss water content dynamics very well, particularly for extended extremely dry periods, because the moss water content can be decoupled from the WTD and soil water content once there is an extended extremely dry period (St-Hilaire et al., 2010; Wu et al., 2011). As St-Hilaire et al. (2010) stated, these aspects need to be included in the MWM once new studies are available that can be used to parameterize these processes in peatland ecosystems.

7 Conclusions

To simulate the sensitivity of C cycling to projected climate change in northern peatlands, the McGill Wetland Model (MWM) was coupled to a land surface climate model (the wetland version of Canadian Land Surface Scheme (CLASS3W)). In this paper, we evaluated this coupled model for a bog (MB) and a fen (DS) by comparing the simulated NEP with the observed NEP from an eddy-covariance tower. CLASS3W-MWM simulates the magnitude, direction, seasonal and interannual variability of CO₂ fluxes reasonably well for both MB and DS, although the summer peak GPP and TER are not reproduced well. The statistical indices for the daily comparison have shown that the performance of CLASS3W-MWM for the simulation of C fluxes in northern peatlands based on the simulated soil and surface climate (r ² = 0.63, 0.42; d* = 0.88, 0.80; RMSE = 0.54, 0.54 for MB and DS, respectively) is comparable to that of the stand-alone MWM based on the observed soil and surface climate (r ² = 0.58, 0.52; d* = 0.80, 0.84; RMSE = 0.21, 0.49 for MB and DS, respectively).

Our sensitivity analysis based on CLASS3-MWM confirms that northern peatlands are thermally and hydrologically conservative ecosystems. There are significant differences between the sensitivity analysis based on CLASS3W-MWM and that based on the stand-alone MWM. Both bogs and fens can change from being C sinks, as they are currently, to being C sources within the range of projected climate changes. This is a result of both increased autotrophic and heterotrophic respiration. But for fens, decreased photosynthesis also plays a role under extreme dry or wet conditions. Changes in air temperature will have more significant impacts on both bogs and fens than changes in precipitation. Fens have a narrower tolerance to environmental changes than bogs. Although both bogs and fens are very vulnerable to drought, which could result from either increased air temperature or decreased precipitation, only changes in precipitation do not cause northern peatlands to become a C source to the atmosphere. Moreover, our sensitivity analysis based on the combination of the individual changes suggests that northern peatlands respond to changes in air temperature, precipitation and CO₂ concentration in a highly non-linear way. The sensitivity of C cycling in northern peatlands with respect to changes in air temperature, precipitation and the concentration of atmospheric CO₂ together is not a simple addition or subtraction of the sensitivity to the individual changes.

Acknowledgements

NTR received support for the development and evaluation of MWM from a Canadian Foundation for Climate and Atmospheric Sciences research network (Canadian Global Coupled Climate Carbon Modelling) and a Natural Sciences and Engineering Research Council of Canada Discovery Grant. MN received support for the micrometeorological instruments and their maintenance at the Degerö Stormyr site from the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning, the Swedish Research Council, and the Kempe Foundation. JW acknowledges the support of the Global Environment and Climate Change Research Centre (GEC3), McGill University (FQRNT Centres Grant) in providing funds to travel to Sweden for part of this research. JW would like extend his thanks to the McGill Recruitment Fellowship for funding his study in his first year at McGill University. JW extends his thanks to the funding support from Grenfell Campus' start-up funding, Principal's Research Fund and Humber River Basin Project fund. We also thank Fluxnet-Canada, now the Canadian Carbon Program (CCP) for providing funds for maintaining the eddy-covariance tower measurements at Mer Bleue. The comments from two anonymous reviewers have greatly improved this manuscript.