Integrating peatlands into the coupled Canadian Land Surface Scheme

10 Peatlands, which contain large carbon stocks that must be accounted for in the global carbon budget, are 11 poorly represented in many earth system models. We integrated peatlands into the coupled Canadian 12 Land Surface Scheme (CLASS) and the Canadian Terrestrial Ecosystem Model (CTEM), which together 13 simulate the fluxes of water, energy and CO2 at the land surface –atmosphere boundary in the family of 14 Canadian Earth System Models (CanESMs). New components and algorithms were added to represent the 15 unique features of peatlands, such as their characteristic ground floor vegetation (mosses), the slow 16 decomposition of carbon in the water-logged soils and the interaction between the water, energy and 17 carbon cycles. This paper presents the modifications introduced into the CLASS-CTEM modelling 18 framework together with site-level evaluations of the model performance for simulated water, energy and 19 carbon fluxes at eight different peatland sites. The simulated daily gross primary production and 20 ecosystem respiration are well correlated with observations, with values of the Pearson correlation 21 coefficient higher than 0.8 and 0.75 respectively. The simulated mean annual net ecosystem production at 22 the eight test sites is 87 g C m yr, which is 22 g C m yr higher than the observed annual mean. The 23 general peatland model compares well with other site-level and regional-level models for peatlands, and 24 is able to represent bogs and fens under a range of climatic and geographical conditions. 25


20
Peatlands represent about 20 % of the global soil carbon (C) pool and have played a critical role in regulating the global climate since the onset of the Holocene . Peatlands have accumulated more than 600 Gt C over the Holocene and serve as a long-term C sink at a rate higher than 5 Gt C per century on average (Yu et al., 2010). Over 90 % of the world's peatlands are located in the Northern Hemi-Siberian Lowlands and the FennoSoviet Lowlands, where gross primary production (GPP) is comparatively low (e.g. Yebra et al., 2015). The inhibited decomposition in waterlogged organic soil persistently sequesters C in peatlands, despite the relatively low primary production.
Peatlands are usually characterized by a ground layer of bryophytes covering 80-5 100 % of the surface (Vitt, 2014). Bryophytes, especially Sphagnum mosses, are nonvascular land plants that are able to effectively capture and store water and nutrients (Turetsky, 2003). Globally, bryophytes and lichens are widely present, especially over tundra, boreal forest floor and desert, and are estimated to account for a net C uptake of 0.34 to 3.3 Gt C yr −1 (Porada et al., 2013), out of 5.0 (±0.9) Gt C yr −1 global net C 10 uptake by land and oceans between 1960 and 2010 (Ballantyne et al., 2012). Peatlands are particularly vulnerable to C loss under climate change. The IPCC Fifth Assessment Report (AR5) projected a large increase of temperature and a risk of lower soil moisture (Christensen et al., 2013;Seneviratne et al., 2010) in the boreal region. Warmer temperatures and drought can both stimulate the decomposition of 15 peat and further enhance climate change through increased CO 2 and CH 4 emissions (Davidson and Janssens et al., 2006;Tarnocai, 2006;Ise et al., 2008;Dorrepaal et al., 2009;Wu and Roulet, 2014). However, the increasing atmospheric CO 2 concentration and temperature may also promote increased primary production and shifts in vegetation ecozones, compensating for the additional C loss from soil respiration (Camill and 20 Clark, 2000;Ward et al., 2013;Wang et al., 2015). Overall, large uncertainties prevail in the future carbon budget of peatlands and its feedback to climate change (McGuire et al., 2009).
Earth system models (ESMs) simulate the global C cycle and feedbacks to climate and are used to make future climate projections. Poor representation of processes re-mate Model (GCM) (Verseghy, 1991;Verseghy et al., 1993), and has been under continuous development since then. It simulates the energy and water balances of the components of the land surface, mainly the temperatures and liquid and frozen water contents of the vegetation, snow and soil for four sub-areas of each grid cell (bare soil, vegetation covered ground, snow covered ground and vegetation over snow). The 10 model has been parameterized for mineral, organic or mixed soil types (Letts et al., 2000). The organic soil parameterization significantly improved the simulations of soil water and energy balances in peatlands and other organic soils (Comer et al., 2000;Bellisario et al., 2010).
CTEM simulates the terrestrial ecosystem C cycle for nine plant functional types 15 (PFTs) and soil through photosynthesis, autotrophic and heterotrophic respiration based on parameterizations developed by Arora (2003) and Arora and Boer (2005). CTEM's treatment of soil moisture and soil carbon pools showed comparatively high correlations with the biome soil pool and turnover time among ESMs (Todd-Brown et al., 2013). These processes determine the flow of carbon in and out of the model's Introduction CTEM directly controls the stomatal activity and the associated stomatal resistance of the PFTs and thus affects the energy and water exchanges at the surface in CLASS. Photosynthesis and leaf respiration are modelled at a time step of 30 min, whereas the rest of the terrestrial ecosystem processes are modelled at a daily time step.
To account for the eco-hydrological and biogeochemical interactions among vegeta- 5 tion, atmosphere and soil in peatlands, the following modifications were made to the coupled CLASS3.6-CTEM2.0 modelling framework: 1. The top soil layer was characterized as a moss layer with a higher heat and hydraulic capacity than a mineral soil layer. The moss layer buffers the exchange of energy and water at the soil surface and regulates the soil temperature and 2. Three peatland vascular PFTs (evergreen shrubs, deciduous shrubs and sedges) as well as mosses were added to the existing 9 CTEM PFTs. These peatlandspecific PFTs are adapted to cold climate and inundated soil with optimized plant structure (shoot/root ratio, rooting depth), growth strategy and metabolic acclima- 15 tions to light, water and temperature.
3. We considered the soil inundation stress on microbial respiration in the litter C pool. The original CTEM assumed that litter respiration was not affected by oxygen deficit as a result of flooding, since litter was always assumed to have access to air. This assumption does not hold for peatlands where high water table positions 20 occur routinely. 4. We separated the soil C balance and heterotrophic respiration (HR) calculations for peatland and non-peatland fractions for each grid cell in the global model. Over the non-peatland fraction, we used the original CTEM approach that aggregated the HR from each PFT weighted by the fractional cover. Over the peatland 25 fraction the soil C pool and decomposition are controlled by the water

Soil layers
The water table depth (WTD) in natural peatlands fluctuates seasonally from above the soil surface to the top of the permanently saturated soil layer, which is often referred to as the boundary between acrotelm and catotelm. The boundary is usually estimated to be 30 cm below the soil surface in wetlands (National Wetland Working Group, 1997), 5 and has been widely used as the bottom of the first soil layer in two-layer soil decomposition models (e.g. Granberg et al., 1999;Yorova et al., 2007;Spahni et al., 2013).
To capture the effect of the fluctuating water table on the transfer of water and energy within the soil, we used a multi-layer configuration rather than the standard three-layer configuration of the soil layers in CLASS. We assigned nine organic soil layers, each 10 10 cm thick, at the top of the soil profile and a 10th soil layer from 90 cm down to the bottom of the organic soil ( Fig. 1). Moss was treated as the top first soil layer and the substrate below the 10th soil layer was considered as bedrock. Mineral soil was not included. 15 The standard configuration of soil layers in CLASS consists of 3 layers with thickness of 0.10, 0.25, and 3.75 m. Organic soil in CLASS was parameterized by Letts et al. (2000) as fibric, hemic and sapric peat in the three soil layers respectively, representing fresh, moderately decomposed and highly decomposed organic matter. Tests of CLASS on peatlands revealed improved performance in the energy simulations for fens and bogs 20 with this organic soil parameterization. However, the model overestimated energy and water fluxes at bog surfaces during dry periods due to the neglect of the moss cover (Comer et al., 2000).

A moss layer as the first soil layer
To take into account the interaction amongst the moss and the soil layers and the overlying atmosphere for energy and water transfer, we added a new soil layer 0.10 m 25 thick above the fibric organic soil to represent living and dead peatland bryophytes, such as Sphagnum mosses and true mosses (Bryopsida). The physical characteris-GMDD 8,2015  tics of mosses differ from those of either the shoots or the roots of vascular plants (Rice et al., 2008). In particular, mosses can hold more than 30 g of water per gram of biomass (Robroek et al., 2009). More than 90 % of the moss leaf volume is occupied by the water-holding hyaline cells (Rice et al., 2008), which retain water even when the water table depth declines to 1-10 m below the surface (Hayward and Clymo, 1982). 5 The parameter values of the moss layer for water and energy properties were derived from a number of recent experiments measuring the hydraulic properties of mosses (Price et al., 2008;Price and Whittington, 2010;McCarter and Price, 2012; Table 1). Living mosses range from 2-3 to over 5 cm in height (Rice et al., 2008) and have lower values of dry bulk density and field capacity than fibric peat (Price et al., 2008).
Compared to fibric peat, the saturated hydraulic conductivity of living moss is higher by orders of magnitude (Price et al., 2008) and the thermal conductivity is more affected by the water content (O'Donnell et al., 2009). To fully account for the effect of mosses, we set the depth of the living moss (d m ) within the top soil (i.e. moss) layer to 3 cm for fens and 4 cm for bogs, and interpolated its water content w m from the water content of 15 the overall layer and the depth of the living moss: where the dry moss biomass (B m ) is converted from moss C (C m ) using the standard conversion factor of 0.46 kg C per kg dry biomass, θ l,1 (m 3 m −3 ) is the liquid water content of the top soil layer, and ρ w is the density of water (1000 kg m −3 ). The maximum 20 and minimum moss water contents were estimated from a number of observed moss water contents (e.g. Flanagan and Williams, 1998;Robroek et al., 2009). In CLASS, evaporation at the soil surface is controlled by a soil evaporation efficiency coefficient β (Verseghy, 2012). This parameter is calculated from the liquid water content and the field capacity of the first soil layer following Lee and Pielke (1992 where w m , w m,max , w m,min are the water content and the maximum and minimum water contents of the living moss in kg water per kg dry moss. 5 Mosses are an important contributor to the primary production and the C sequestration in peatlands, owing to the low decomposability of the moss tissue. Sphagnum in peatlands grows at 20-1600 g biomass m −2 yr −1 and accounts for about 50 % of the total peat volume (Turetsky, 2003). We have modified CTEM to include a moss C pool and moss litter pool along with the related C fluxes, i.e. photosynthesis, autotrophic 10 respiration, heterotrophic respiration and humification. The net photosynthesis of moss (G m ) is calculated from the gross photosynthesis (G 0,m ) and dark respiration (R d,m ):

Primary production of mosses
The moss photosynthesis and dark respiration are calculated using the Farquhar (1985) biochemical approach following the MWM (St-Hilaire et al., 2010) and CTEM 15 (Melton and Arora, 2015), with modifications for integration with CLASS-CTEM and moss phenology. The leaf-level gross photosynthesis rate G 0,m (µmol CO 2 m −2 s −1 ) is obtained as the minimum of the transportation limited photosynthesis rate (J s ) and the first root of the quadratic solution of the light-limited rate (J e ) and the Rubisco limited rate (J c ). A logistic factor (ς) is added with values 0 or 1 to introduce a seasonal con-20 trol of moss photosynthesis. In the MWM, spring photosynthesis starts when the snow depth is below 0.05 m and the soil temperature at 5 cm depth goes above 0.5 • C . Since in our case CLASS sets the minimum depth for melting, discontinuous snow to 0.10 m, this limits the spring photosynthesis to starting only once the snow 10097 Introduction The dark respiration in mosses (R d,m ) is calculated as a function of the base dark respiration rate (R d,m,0 ) which has a value of 1.1 µmol m −2 s −1 (Adkinson and Humphreys, 2011) scaled by the moss moisture (f m,rd ) and soil temperature functions (f T ,rd ). The

5
MWM models the relation between water content in mosses and dark respiration with optimal water content at 5.8 g water per g dry weight, following the approach in Frolking et al. (1996). We modified the relation for water content above the optimal water content, based on a recent discovery of a weak linear positive relation between the dark respiration rate and the water content above the optimal water content during the late Photosynthetic photon flux density (PPFD) is measured by the photosynthetically active 15 radiation (PAR), which is defined as the solar radiation between 0.4 to 0.7 µmol that can be used by plants via photosynthesis. In the coupled CLASS-CTEM system, the PAR received by the moss (PAR m , unit µmol protons m −2 s −1 ) is converted from the visible short-wave radiation reaching the ground (K * g , unit: W m −2 ) in CLASS by a factor of 4.6 µmol m −2 s −1 (McCree, 1972). K * g is a function of the incoming shortwave radiation (K ↓, unit: W m −2 ), the surface albedo (α g ), and the canopy transmissivity (τ c ): The energy uptake by the moss layer is thus a function of the total incoming short-wave radiation, the aggregated leaf area index (LAI) of the PFTs present, the snow depth, the fractional vegetation cover and the soil water content (Verseghy, 2012). In peatland C 5 models that do not consider vegetation dynamics, the transmissivity of the vegetation canopy is usually assumed to be constant (e.g. St-Hilaire et al., 2010). Compared with such models, CLASS enables a better representation of light incident on the moss surface since it includes partitioning of direct/diffuse and visible/near-IR radiation, PFTspecific transmissivities, and time-varying LAI and fractional PFT coverages (Verseghy et al., 2012).  et al., 2003) and they share similar responses to climate in ESMs (e.g. Bonan et al., 2002). Table 2 lists the key parameters for the peatland PFTs used in this model. (The photosynthesis and autotrophic respiration of vascular PFTs are modeled the same as the original CTEM.) 5 Over the non-peatland fraction, heterotrophic respiration (HR) is calculated as the sum of the respiration from litter and soil carbon pools as in the original version of CTEM (Arora, 2003). The soil C pool over the non-peatland areas is assumed to be exponentially distributed with depth (Arora, 2003). In peatlands a large amount of humic soil is generally located in the permanently saturated zone and the bulk density increases 10 with soil depth (Loisel and Garneau, 2010). Thus the assumption of exponentially decreasing distribution of C content with increasing soil depth is not valid in peatlands. We used a quadratic equation to calculate the distribution of soil C content over depth based on an empirically determined bulk density profile (Frolking et al., 2001). HR over the peatland fraction of a grid cell is modelled using a two-pool approach 15 with a flexible boundary between the pools that depends on the depth of the water table:

Heterotrophic respiration
where o and a denote the oxic and anoxic portions of the soil C pool, respectively. The respiration rate R (unit: µmol C m −2 s −1 ) is obtained from the respiration rate constant k 20 (µmol C kg C −1 s −1 ), the temperature functions f T , the soil C mass C SOM (kg) and a scaling factor f anoxic which is set to 0.025 (Frolking et al., 2010) representing the inhibition of microbial respiration under anoxic conditions. Q 10 is calculated using a hyperbolic tan function of the soil temperatures (T s ) of the oxic and anoxic zones (Melton and Arora, 2015), which are in turn functions of water  et al. (2001), and parameterized differently for fens and bogs (Table 3): 10 where 0.487 is a parameter that converts from soil mass to soil C content. As only organic soil is considered in peatlands, the peat soil C is updated from the humification (C hum kg C m −2 day −1 ) and soil respiration from the oxic (R o in kg C m −2 day −1 ) and GMDD 8,2015 Integrating peatlands into the coupled Canadian Land Surface Scheme anoxic (R a in kg C m −2 day −1 ) components during the time step: At the end of each time step, the peat depth (i.e. the depth of the organic soil) d p is updated from the updated peat C mass (C SOM in kg) by solving the quadratic equation: The water table depth z wt is deduced by searching for a soil layer below which the soil is saturated and above which the soil moisture is at or below the retention capacity with respect to gravitational drainage. Within this soil layer j , z wt is calculated as: where ∆z is the thickness of soil layer (unit: m), θ l and θ i are the liquid and frozen water 10 contents (unit, m 3 m −3 ), θ ret and θ p are the water retention capacity and the porosity, and z b (unit: m) is the bottom depth of the soil layer.

Site locations
The model was applied at eight peatland sites to assess its performance in simulat-  (Table 4).

Model initialization and spin up
For each site, the FLUXNET database was used to assign values to background variables such as latitude, longitude, peat depth, areal coverage of the three peatland PFTs, and their roughness lengths, visible and near-infrared albedos and canopy mass.
Other CLASS-and CTEM-related vegetation parameters were assigned their standard values, as listed in Table 2. The parameter values for evergreen shrubs, deciduous shrubs and sedge mostly reflected those used for evergreen needleleaf trees, deciduous needleleaf trees and C3 grasses in CTEM, respectively. Exceptions were made for some parameters that determine the length or shape and turnover of the stem and root of the PFT and its tolerance to coldness and dryness (Table 2). 20 Model pools were spun up from initial conditions by repeatedly cycling through the inputs for approximately 100 years until the annual mean C pools in vegetation in consecutive years differed by less than 5 %. The initial soil C mass was calculated from the observation-based estimations of peat depth based on an empirically obtained relation between the soil depth and soil mass (Eq. 15). 8,2015 Integrating peatlands into the coupled Canadian Land Surface Scheme

Observational data sets
The model was forced with half-hourly measured meteorological data: downwelling shortwave radiation, downwelling longwave radiation, precipitation, atmospheric pressure, air temperature (T a ), specific humidity, and wind speed. The measurement heights for the latter three were obtained from the FLUXNET metadata. Datasets ranged in 5 length from 2 to 9 years. The parameters used for model evaluation include water table depth (z WT ), snow depth, soil temperature (T s ), latent heat flux (QE), sensible heat flux (QH), GPP, ER and NEP. Energy and C fluxes were measured every 30 min using the eddy-covariance (EC) technique. The required downwelling longwave radiation (LW ↓) was available only at MB-Bog, AB-Fen, SE-Deg and FI-Lom. For the remaining 4 sites, 10 LW ↓ was estimated following the methods of Crawford and Duchon (1998): where σ is the Stefan-Boltzmann constant and c f is the cloud fraction term ranging between 0 and 1. c f is estimated as the ratio between the incoming shortwave radiation and the clear-sky solar radiation, which in turn is a function of the locational character 15 of the site, i.e. latitude, longitude, altitude and time zone. ε c is the clear sky emissivity and is estimated from the vapor pressure (e 0 ) following Ångström (1918): Water

Evaluation methods
The model was evaluated against observation-based sensible and latent heat fluxes at the soil surface, soil water content, water table and snow depth, soil temperature at various depths and the daily, monthly and annual C fluxes (GPP, ER, NEP). The root mean square error (RMSE) and linear regression coefficient (r 2 ) were primarily used 5 for evaluation. Statistical analyses were conducted using the free software package R version 3.1.1 (R Core Team, 2014).
Since the ultimate goal is to apply the model globally in an ESM, further experiments were done to investigate the importance of modelling fens and bogs separately. In the version of the model described above, bogs and fens are distinguished primarily through the parameterization of the control of water table depth on soil decomposition (Table 3). Also, the depth of the living moss (d m ) is set to 4.0 cm for bogs and 3.0 cm for fens. In a first test, the parameters for soil decomposition (Table 3) for bogs were used for the fen sites and those for the fens were used for the bog sites. In a second test, the living moss layer was set to a set to a single fixed value of 3.5 cm for both bogs 15 and fens. The resulting differences in the surface fluxes and the soil temperatures were then evaluated.  (Table 1). Correcting the modeled WTD by 0.25 m led to a high agreement with the observed WTD in MB-Bog (Fig. 3). For AB-Fen, the model overestimated the inter-annual fluctuation and did not reproduce the trend of increasing WTD seen in the observations, which was likely associated with the change in vegetation 15 cover. It has been observed that the AB-Fen site is currently changing from a rich fen to a poor fen and is now in a phase of rapid tree establishment and increase in LAI and NEP . The model reproduced the annual variation of snow depth quite well for the bog and fen sites where observations were available (Fig. 4). The errors for the MB bog may be 20 associated with uncertainties in the observed data stemming from the combination of a continuous record from one spot with sporadic snow depth data from other locations on the bog surface .

Energy budget terms
The model performed similarly well on the daily latent heat (QE) and sensible heat 25 (QH) fluxes for multi-year simulations ( that the simulation periods for our sites ranged from 2 to 9 years compared to the 1 year simulation with LPJ-WHy, and that we included eight sites in our evaluation compared with two peatland sites for LPJ-WHy. Our model was able to capture the seasonal variation in soil temperature at different depths down to the bedrock. Figure 6 compares the modeled soil temperatures against the observations at 5, 40, 80, and

Carbon fluxes
Examination of the modelled GPP, ER and NEP demonstrates that the model is capable of capturing seasonal dynamics and climate-driven events consistently in various types of peatlands . For example, the RU-Fyo bog experienced a period of low GPP due to an abrupt decrease of air temperature in the early fall of 2010, which 5 was well reproduced by the model. The RMSE ( Table 6)   observed NEP varied widely, ranging from −1.8 to 2.2 g C m −2 day −1 and from −3.9 to 4.8 g C m −2 day −1 respectively. Model errors for the extreme values at these two sites may have contributed to their low r 2 values. NEP was overestimated at the beginning and the end of the growing season for FI-Lom due to the overestimation of GPP for that period as discussed above. These results may be compared to an evaluation of 5 the MWM using the SE-Deg dataset that was conducted by Wu et al. (2013). For daily NEP they obtained an RMSE of 0.49, similar to ours, but a higher r 2 of 0.52. It should be noted that the MWM was driven by observed WTD and soil temperature, while in our simulations these were allowed to evolve freely, so our reasonable result is encouraging. 10 The simulated accumulated monthly NEP from March to November agreed well with the observations in the four bogs and four fens. The outliers for bogs were the overestimations in MB-Bog in October and November due to the underestimation of GPP (Fig. 7). The NEP in RU-Fyo in one August was underestimated owing to the underestimated GPP, which in turn was a result of the underestimated LAI and rooting depth 15 temperature in the summer. Figure 10, showing plots of NEP averaged for each month of the year at each site, demonstrates on the whole larger scatter for the bogs than the fens, with the scatter increasing through the summer and fall. The overall value of r 2 was 0.59 for bogs and 0.58 for fens; both values are higher than or similar to those obtained in evaluations of other peatland C models. For example, the r 2 value of the

Annual carbon budget
The simulated mean annual NEP values with their standard deviations generally fall within the range of the standard deviations of the observations (Fig. 11), between 9 g C m −2 yr −1 in the rich fen (FI-Lom) and 73 g C m −2 yr −1 in the productive bog (RU-Fyo) ( Table 7). The only site with large bias in annual NEP was AB-Fen. Observation-5 based estimations of NEP in this fen were extremely high, totalling 176 g C from May to October, in comparison with other sites (Syed et al., 2006). This treed fen had a high peat density and LAI and large variation in the WTD, which, accompanied by high spring temperatures, resulted in high ecosystem photosynthesis capacity and production (Adkinson et al., 2010). Considering nutrient factors and the site-specific peat den-10 sity could potentially capture the large NEP at this site. The observed annual NEP for the eight sites varied greatly overall, between −17 and 187 g C m −2 yr −1 , while the simulated NEP showed slightly less variation, ranging from 13 to 157 g C m −2 yr −1 . The simulated mean annual NEP across the sites was 87 g C m −2 yr −1 and was 22 g C m −2 yr  (Dinsmore et al., 2010) while the dry MB-Bog was estimated to be a source of 13.8 g C m −2 yr −1 (Roulet et al., 2007). The modeled NEP bias tended towards underestimation for the treed fen (AB-fen) and the productive ombrotrophic bog (MB-Bog), and towards overestimation for the remaining sites (Fig. 11).

25
The model errors in GPP were smaller than the standard deviation of the observations, except for the atypical sites (AB-Fen, RU-Fyo) and the sites that had only a few years of data (FI-Lom, SE-Faj) ( the error bars except for in the RU-Fyo bog, for which a thin peat depth of 1 m was used to initialize the simulation ( Table 4). The simulated WTD was consistently shallower in the summer than the observations (Fig. 3), which slowed down the soil respiration in the model and contributed to the discrepancies in ER. The observed WTD showed an abrupt decrease in the summer of 2010 without pulses of large ER being observed 5 during that period (Fig. 8), indicating uncertainties in the WTD observations. Another reason for the errors in ER was the underestimation in soil T . For example, the simulated soil T at 5 cm depth was higher in the summers with RMSE of 4.6 • C in RU-Fyo (Table 5). The site is particularly shallow and homogeneous, thus the standardized living moss layer of 4 cm for bogs was probably too large, leading to an overestimation of the thermal insulation effect from the moss layers and hence less seasonal variation in soil temperature and ER. An overview of the model's performance is illustrated via a Taylor diagram (Fig. 12). This demonstrates the model's skill in simulating the hydrological and thermal dynamics and C fluxes in different types of peatlands across a variety of climatic and geographical 15 settings. The outliers are the vegetated treed fen (AB-Fen), the maritime blanket bog UK-Amo and the extremely shallow peatland RU-Fyo. The model simulations consistently agreed quite well with the observations except at these sites for some evaluated parameters. The Pearson r was above 0.90 for the soil temperature at 5 cm and above 0.50 and 0.60 for the sensible and latent heat fluxes, except for those at UK-Amo. The 20 soil water content at the surface soil layer, corresponding to the available measurements of soil temperature at 5 cm depth, was about 0.6 with RMSE between 0.20 and 0.28 m 3 m −3 . The modeled daily GPP and ER were highly correlated with the observations, with Pearson r values between 0.80 and 0.95 for GPP, and between 0.75 and 0.96 for ER. The simulated daily NEP accumulated the errors in GPP and ER and was somewhat less well correlated with the observations, with Pearson r values between 0.4 and 0.72. 8,2015 Integrating peatlands into the coupled Canadian Land Surface Scheme

The necessity of distinguishing fens and bogs
The original version of our peatland model (referred to as "CONTROL" hereafter) as described above distinguishes bogs and fens through the controls of water table depth on soil decomposition and the depth of the living moss. The parameters for the water table depth regulation of soil decomposition were derived from the empirical relations 5 in the MWM (Eqs. 13 and 14). Our first test, "K-SWAP", involved swapping the values of the decomposition parameters (Table 3) between the bog and fen sites. As shown in Fig. 13, the differences between the test and control runs are minimal. The relative differences in the simulated values of the fluxes and temperatures between K-SWAP and CONTROL ranged from −1.6 to +5.1 % for RMSE and from −23 to +6 % for r 2 . The 10 relative differences in RMSE and r 2 for GPP, QH, QE and T s5 were smaller than ±1 %.
The largest differences in r 2 between K-SWAP and CONTROL were for NEP at SE-Faj and UK-Amo, which had significantly lower r 2 values than the other sites. The results of K-SWAP indicate that parameterizing fens and bogs differently for the regulation of water table depth on soil decomposition makes little difference in the simulation. 15 The second test, "D-MOSS", retained the settings in K-SWAP and changed additionally the depth of the living moss in both bogs and fens to 3.5 cm. The RMSE and r 2 of D-MOSS show site-specific differences compared to CONTROL (Fig. 13). The relative differences between D-MOSS and CONTROL in RMSE and r 2 were in the range of −5 to +7 % and −15 to +13 %, respectively. The mean differences for all sites and all 20 evaluated variables were less than 5 % for both RMSE and r 2 . For GPP, ER and the soil temperature at 5 cm depth, the r 2 in D-MOSS was similar to that of CONTROL. For QE, the r 2 in D-MOSS was higher than the control for all the fens and one unusual bog (UK-Amo), but not for the other three bogs. Compared to CONTROL, the r 2 of NEP was higher in D-MOSS for five sites by up to 7 % and less than 2 % lower 25 in the other sites, except for UK-Amo where r 2 was also low in CONTROL. This test indicates that the depth of living moss is important for modeling energy exchange and net C exchange at the ground surface. However, the depth of living moss seems to be more site-specific than the bog vs. fen difference. The two tests above suggest when applying the model, it is not necessary to distinguish between fens and bogs, in contrast for example to the MWM and its soil decomposition component the Peatland Decomposition Model (PDM) (Frolking et al., 2001), 5 which were developed for the detailed modelling of specific sites. Therefore, when the present model is implemented within CLASS-CTEM on regional and global scales, one general type of peatland may be simulated with no differentiation between bogs and fens. This will considerably simplify the global implementation of the model, since global datasets mapping the locations of fens vs. bogs are not available. 10

Conclusions
We have presented here an extension of the CLASS-CTEM model, enabling it to simulate the water, energy and C cycles of peatlands. The model simulations of the daily C fluxes are of comparable accuracy to those performed by other models that were developed for a particular site or an area, for example the Finland regional peatland model 15 (Gong et al., 2013) for the FI-Lom site and the MWM for the MB-Bog and SE-Deg sites . Compared with models that simulate global peatland C fluxes such as LPJ-WHy (Wania et al., 2009a, b) and CLIMBER2-LPJ (Kleinen et al., 2012), our model performs well and covers the ranges in the observations (Yu et al., 2010). The variations in climatic conditions and in the C stocks contained by peatlands in nature 20 are difficult to capture completely by the general peatland model here. The model errors were larger for sites with unusual soil properties or vegetation cover. Long-term decline of water table depth can also shift the vegetation in peatlands from mosses and grasses to shrubs and trees Munir et al., 2014;Talbot et al., 2010). Taking into account such effects could improve the performance of the 25 model (Sulman et al., 2012). Also, other forms of C besides CO 2 , such as methane (CH 4 ) and dissolved organic C, are as yet missing from the C budget in the model and need to be included in order to fully simulate the net C budget of peatland ecosystems. At the moment, approaches to modelling CH 4 emissions from peatlands or wetlands diverge widely and further work is needed in areas such as more accurate land surface classification, more realistic emissions from non-inundated wetlands (where water table depth regulates the emissions) and peat soils from high latitudes (Bohn et al., . This study has tested the model's performance on northern peatlands only; further tests are needed to validate the model on the remaining 10 % of peatlands (Yu et al., 2011) that are located in the tropical region and Southern Hemisphere. The coupled CLASS-CTEM models serve as the land surface component for the family of Canadian Earth System Models (CanESMs). Despite some limitations in simulating unusual peatlands, the extended version that we have presented here shows an overall good skill in simulating the water and energy dynamics and the daily and annual C fluxes in peatlands. Contrary to models designed for specific sites such as the MWM, the peatland model presented here need not distinguish between bogs and fens, which constitutes a distinct advantage for application in an ESM at the global 15 scale.