A new model of Holocene peatland net primary production , decomposition , water balance , and peat accumulation

Introduction Conclusions References


Introduction
Northern peatlands contain somewhere between 250 and 600 Pg C (Gorham, 1991;Turunen et al., 2002;McGuire et al., 2009). At present, multi-year field studies indicate that northern peatlands take up CO 2 at a rate of 40-80 g C m −2 yr −1 and lose 10-20 g DOC m −2 yr −1 to downstream aquatic ecosystems Nilsson et al., 2008;Koehler et al., 2010), while northern peatlands in aggregate 15-40 Tg CH 4 yr −1 (Makiloff-Fletcher et al., 2004). Peatlands have played an important role in the greenhouse gas composition in the atmosphere for most of the Holocene, leading to an estimated net cooling of ∼0.5 W m −2 (Frolking and . Paleo-ecological and biogeochemical studies have conjectured the role of peatlands in the Holocene dynamics of atmospheric CO 2 and CH 4 (MacDonald et al., 2006). Despite the possible global influence in climate over the Holocene, peatlands are just beginning to be incorporated into Holocene climate assessments. A few studies have hard-coded peatland carbon sinks into Holocene climate simulations (e.g. Wang et al., 2009;Kleinen et al., 2010) but in these assessments there are no (Wang et al., 2009) or only one-way (Kleinen et al., 2010) interactions between climate and peatland carbon accumulation. There is therefore a need to develop relatively simple yet realistic peatland growth models where the rates of ecosystem processes are a function of climate and the inherent autogenic properties of peatlands (Belyea and Baird, 2006;Dise, 2009). In the present study we develop a new model, with the intention that, with further development (see Discussion in Sect. 5 below) it could be eventually coupled to Holocene climate or earth system models to investigate the climate-carbonmethane interactions over decades to millennia.
There have been some attempts to simulate peatland development, particularly the growth of peat bogs. Most peatland modeling stems from the seminal work of Clymo (1984), and the model developed in this paper has its roots in that work. Our objective is to develop a simple, one-dimensional model of peat accumulation that explicitly includes the feedbacks among hydrology, plant communities, and peat properties, begin to evaluate how these feedbacks affect peatland development history, and test if this model structure is able to simulate basic patterns of northern peatlands peat accumulation over millennia, including varying accumulation rates, and vegetation and fen-bog transitions. Longer term goals include (1) providing insight into climate domain(s) in which northern peatlands are persistent accumulating systems, (2) making projections of impacts of future climate change on northern peatland carbon and water dynamics, and (3) contributing to a methodology for incorporating peatlands as unique land surface types into earth system models ).
2 Background-modeling peat processes Clymo (1984) developed the basic concept that underlies much of our ideas of long-term peat accumulation and the interpretation of peat core data as a record of carbon accumulation (there is a much longer history of peat core analysis as a record of paleoecology and paleoclimate). In the Clymo model, a peat profile is viewed as two-layered -an overlying acrotelm above the long-term mean low (summer) water table (Ingram, 1982) which receives plant litter (mass) inputs, and an underlying, persistently-saturated catotelm layer which receives mass inputs only from the overlying acrotelm. Each layer loses mass through decomposition. Under relatively stable climate conditions, in this model the acrotelm relatively quickly reaches a steady state, and effectively "floats" above a thickening catotelm, and so the Clymo model focuses more explicitly on the longer-term dynamics of the catotelm, which is the much larger peat/carbon reservoir in most peatlands. The change in catotelm mass (M c ) is modeled as where p c is the annual mass transfer across the acrotelmcatotelm boundary, nominally ∼10% of annual net primary productivity or NPP (Clymo 1984), and k c is the annual fractional mass loss rate due to catotelm decomposition (α in Clymo's notation). Integrating Eq. (1) gives that has an equilibrium mass for large t of M * c = p c k c . Inverting this equation to to generate a relationship between cumulative mass per unit area and age, or using a constant value for peat bulk density to convert mass to depth, this model generates a concave depth-age profile. The model does not include any representation of vegetation, and treats water dynamics only implicitly, as the water table must rise each year at the same rate that the catotelm grows to keep a constant water table depth and acrotelm thickness (Belyea and Baird, 2006). In most applications, Eq.
(2) is fit to field data -dated peat core profiles of age (t) vs. accumulated mass (M c ) -with fitting parameters p c and k c , neither of which are directly observable. Further developments with this model have included (i) allowing k c to decrease with increasing mass loss, i.e., as M c (t)/(p c · t) decreases (Clymo 1992, Clymo et al., 1998, which changes the equilibrium mass, and (ii) allowing p c to vary with time, which can generate a convex, linear, or concave age-depth profile (Yu et al., 2003;Belyea and Malmer, 2004). These variations on the model in Eq. (1) also simulate only mass accumulation in the catotelm, and ignore or only implicitly include vegetation and hydrological dynamics. Operationally, all of these models are used to fit an equation to profile data to find parameter values; often several models can generate statistically equivalent fits to observed data (Clymo, 1992;Clymo et al., 1998). Frolking et al. (2001) developed a simple bookkeeping model that tracked annual litter cohorts of two plant functional types (vascular and non-vascular) as they were buried and, through time, "moved" down the peat profile. This model showed that a plausible peat profile could be generated using observable parameters -vegetation net primary productivity (NPP; e.g., Wieder 2006 and references therein), aboveground-belowground partitioning of NPP for the vascular vegetation (e.g., Rydin and Jeglum, 2006;Murphy, 2009), initial (surface) litter decomposition rates from litter bags studies (e.g., Moore et al., 2007) for each plant type, and observed bulk densities. Belowground productivity was added to litter cohorts within the prescribed rooting zone. S. Frolking et al.: A new model of Holocene peatland net primary production 3 Litter/peat decomposition rates declined linearly with mass lost as where m o is the initial (fresh) litter mass, m is the mass remaining in the litter/peat cohort, and k o is the initial (surface) decomposition rate. This is one of the functions for declining decomposition rates evaluated by Clymo (1992) and Clymo et al. (1998). When the cohort was buried below the prescribed, constant water table depth, decomposition rates were reduced by about an order of magnitude to represent slower anaerobic decomposition -this transition was prescribed to be more rapid for bog than for fen simulations. In this model, assuming a constant water table and NPP, about 90% of the initial mass was lost during the several decades required to "traverse" the acrotelm, and decomposition rates deep in the catotelm were two to three orders of magnitude lower than the initial or surface rates. Litter/peat cohort composition (i.e., the fraction of vascular vs. non-vascular mass remaining) varied down the profile and depended on k o values, water table depth, and above-belowground partitioning of vascular NPP. Bauer (2004) extended this model to include additional plant functional types, litter chemical fractions (soluble compounds, holocellulose, and lignin) with different decay rates, and externally imposed dynamics in water table depth and pH that changed total NPP and the relative productivity of different plant functional types. Bauer's results showed that the interactions of vegetation properties (e.g., litter quality) and environment conditions (e.g., water table depth and pH) can have complex effects on peat accumulation rates that may be difficult to interpret from peat core data. Carbon and water cycling in peatlands are tightly coupled. Peat water content (water table depth and degree of saturation above the water table) affects productivity and decomposition, and the carbon balance and vegetation composition affect peat hydraulic properties and water. Hilbert et al. (2000) modeled peat accumulation at an annual time step as a dynamic interaction between accumulating peat and accumulating water, using two coupled equations where h PD is the total peat height (catotelm plus acrotelm), G is the annual plant litter input, k a and k c are the fractional height annual loss rates due to decomposition in the acrotelm and catotelm (equivalent to a mass loss rate for a given bulk density), respectively, z WT is the water table depth (positive down from the top of the acrotelm, θ max is the water content at saturation, and dW/dt is the annual water balance (precipitation plus run-on minus evapotranspiration minus run-off; dW/dt = P + R on − ET − R off ). By making G and ET simple non-linear functions of z WT , and R off a linear function of h PD , the coupled equations generated dynamical behavior and equilibrium states with deep peat and relatively deep water tables or shallower peats and shallower water tables. This coupling of the carbon and water cycles in a peatland provides a simple representation of a dynamical link that is considered an important factor in peatland autogenesis. These feedbacks can influence the vegetation composition of the peatland which in Eq. (5) was represented by a parabolic function G(z WT ), can lead to a dynamic and variable acrotelm thickness, can affect the composition and physical properties of the litter/peat (pore size distribution, bulk density, hydraulic conductivity, and water retention), and through this can exert important controls on long-term peat accumulation (Belyea and Baird, 2006). The dominant control explored by Hilbert et al. (2000) was the site water balance (precipitation minus maximum evapotranspiration). The model included a representation of the acrotelm, a dynamic annual water balance and aggregate vegetation productivity. Hilbert et al. (2000) did not compare model results to field data, but used it as a theoretical tool to study peatlands as a dynamical system by focusing on the non-linear interactions among peat production, decomposition and hydrology. Recently, the dynamical approach to modelling has been extended to look at the development of patterning on peatlands through the inclusion of the regulating feedbacks of hydrology, nutrients, and plant community structure (e.g., Glaser et al., 1992;Rietkerk et al., 2004;Eppinga et al., 2008Eppinga et al., , 2009. The links between vegetation productivity, litter/peat decomposition, peat properties, and hydrology generate a set of interconnected feedbacks that control the dynamics of the system (Fig. 1a) and make prediction a challenge. Here we develop a simple model that represents some of the basic ideas that have emerged as to how peatlands functionvegetation productivity and its rate controls, aboveground litter inputs as well-stratified layers, root litter inputs into the upper peat profile, litter/peat decomposition and its rate controls, internal and external controls on the water balance and water table depth and vegetation and decomposition responses to this -and evaluate the consequences of the model formulation on long-term (millennial) peatland development and carbon accumulation.

Model description and simulations
The Holocene Peat Model (HPM) combines the peat decomposition model (Frolking et al., 2001) and the dynamic peat accumulation model (Hilbert et al., 2000). HPM is a onedimensional (vertical), annual timestep model that simulates vegetation litter production, litter/peat decomposition and net peat accumulation, peat physical and hydrological properties, and the annual water balance, water table location, and unsaturated zone water content (Fig. 1b) 1. (a) Connections between peatland carbon and water cycles and some links to the climate system; (b) schematic of HPM. above), but HPM simulates peat carbon in addition to peat height, and includes refinements to the productivity and decomposition functions and peat hydrological properties -i.e., each term on the right-hand side of Eq. (5).
The vegetation sub-model in HPM is based on the assumption that the controls that the peat environment exerts on peatland vegetation composition and productivity (NPP per unit ground area) can be adequately described by two factors: (1) annual water table depth (z WT , measured down from the surface), which can vary rapidly (annually or every timestep in HPM; e.g. Roulet et al., 2007) and have or not have a longerterm trend, and (2) total peat depth (h PD ) as a proxy for ombrotrophy, i.e., access to mineral nutrients, or buffering capacity against organic acidity generated by decomposition, which can only vary slowly (except for major disturbances like fire or harvest), and will generally have an increasing trend over time. In this first development of HPM, we are assuming no groundwater exchange so the chemical influence through depth is by diffusion only, but in principle advection could be added.  Table 1 for the parameter values for all 12 PFTs. Note that for vascular PFTs, the water table depth used is the mean of the current and past 10 years, while for bryophytes the current year water table depth is used. Color bar is linear in all panels, but represents a relative scale from low (dark blue) to high (dark red). NPP for all PFTs is scaled by a single value so that the maximum NPP in panel (d) equals a sitespecific prescribed value. Each year the developing peatland has particular values for z WT and h PD , which determine NPP for each PFT. All NPP is deposited as litter (i.e., no live biomass accumulation), and vascular PFT NPP is partitioned into above and belowground fractions (see Table 1 and text for details).

Carbon balance equations
In HPM, peatland vegetation is aggregated (or disaggregated, depending on your world-view) into 12 plant functional types (PFTs). The PFTs are distinguished by their productivity characteristics (relative maximum NPP, optimal peat depth and water table depth, sensitivity to non-optimal peat depth and water table depth), their rooting characteristics (belowground fraction of NPP, root density profile); and their litter tissue quality (Table 1). Seven PFTs represent the vascular plants (two of these are woody shrubs -trees will be added in a later version of HPM), and five PFTs represent the bryophytes. Vegetation composition is determined by the relative productivities of each PFT. Seedling establishment/PFT recruitment is not modeled; all PFTs are always present, though with near-zero productivity in non-optimal conditions (Fig. 2). Annual NPP is modeled as two-dimensional, asymmetric Gaussian functions (Fig. 2) for each PFT.  Fig. 3b) is scaled to a maximum of 3 kg m −2 y −1 . * * k 0 are the tissue initial decomposition rates, determined by k 0 = k(1+3k), where k is first order exponential decay rate, i.e., m(t) = m 0 e −kt (e.g., Moore et al. 2007), and k 0 is the value that gives similar m(t) at t ∼5 years (i.e., near the end of reliable values from a litter bag field decomposition study) when m(t) = m 0 /(1 + k 0 t) (e.g. Frolking et al., 2002), the formulation used in HPM., E.g., if k = 0.15 y −1 , k 0 = 0.22 y −1 , e −5k = 0.47; 1/(1+5k 0 ) = 0.48.
where z opt WT is the optimal annual water table depth and h opt PD is the optimal peat depth for PFT i productivity, and σ ± . Annual litter inputs are equal to annual NPP (i.e., PFT biomass accumulation is not modeled) with aboveground NPP forming a new annual litter cohort on the surface of the peat, and belowground NPP added to sub-surface litter cohorts proportional to root density profile and sub-surface litter cohort thickness (Frolking et al., 2001;Bauer, 2004). Root density profiles follow Bauer (2004); sedge root profiles having an exponential decay and 80% of root litter input occurring in the top 0.3 m of litter/peat, and all other vascular plants have a constant root density profile to the maximum of the annual water table depth or 0.2 m. Peat depth/nutrient chemistry impact on NPP is modeled the same way as water table depth, though in reality it may have more effect on plant abundance than directly on NPP, as an ombrotrophic PFT might do fine in minerotrophic conditions, but would be out-competed by minerotrophic species, so would have low abundance and therefore low NPP per unit ground area.
Each annual litter/peat cohort is an assemblage of 12 litter types, each with different initial masses (aboveground NPP), different root litter additions, and different initial decay rates (Table 1). Litter/peat decomposes via a simple non-linear model, following Frolking et al. (2001), where M (= i m i ), the annual litter cohort mass remaining, is the sum of the individual PFT litter masses, l i is the annual litter input for PFT i , which in the field will include living plant shoots that are overgrown or engulfed by growing moss (e.g., Forrest and Smith. 1975), k i is the initial litter decomposition rate for PFT i , m i,o is the total fresh litter input for PFT i into the cohort, and f n is a scalar multiplier for the effect on decomposition rate of litter/peat water content (W i , as degree of saturation, n = 1,0.3 ≤ f 1 ≤ 1, see Eq. 8 below) or depth below the water table (n = 2, 0.001 ≤ f 2 ≤ 0.3, 0.3 m (9) scale length for anaerobic effect on decomposition rate f min 0.001 -(9) decomposition rate reduction factor at 'full & persistent' anoxia ρ min 50 kg m −3 (10) minimum peat bulk density ρ 70 kg m −3 (10) maximum potential increase in peat bulk density c 3 0.2 -(10) µ value at which bulk density has increased halfway from min to max (Fig. 4) c 4 0.05 -(10) controls steepness of bulk density transition curve ( Fig. 4) 0.5 -(11) factor for minimum annual ET at low water table (Lafleur et al., 2005b)  (18) organic matter particle bulk density (Van Wijk, 1966) see Eq. 9 below). Above the water table, in unsaturated litter/peat, this is given by where W opt is the optimum water content for decomposition and c 1 is a parameter set so that f 1 (1) = 0.3 (non-vegetation parameter values are listed in Table 2). In the saturated zone below the water table, the multiplier, f 2 , drops exponentially with depth below the water table (ẑ = z − z WT ;ẑ > 0), representing effects on decomposition rates of water residence time (end-product accumulation, extreme anoxia; Beer and Blodau, 2007;Beer et al., 2008), as where c 2 is the scale length for this decline in decomposition rate below the water table, and f min is the minimum decomposition multiplier (Frolking et al., 2001). HPM does not consider interactions between litter types, so the long-term decomposition rate of litter from a particular PFT depends on that PFT's m i /m i,o , but not on values for any other PFTs (Hoorens et al., 2010). Loss of litter/peat structural integrity due to decomposition, coupled with accumulating mass above a decay cohort, leads to a collapse/compression of the peat (Clymo, 1991). Peat bulk density can increase by 100% or more, pore size distribution changes markedly, causing peat hydrological properties to change as well -well-decomposed or humified peat generally has orders of magnitude lower values of hydraulic conductivity than fresh litter/peat, and a tendency to hold much more water under the same tension (e.g.,

Fig. 3.
Peat cohort bulk density (ρ; Eq. 10) as a function of degree of decomposition (as fraction of initial mass inputs remaining, m/m 0 ). This curve uses the parameterization in Table 2, with a minimum bulk density of 50 kg m −3 and a maximum bulk density of 120 kg m −3 .
Walmsey, 1977; Paavilainen and Päivänen, 1995;Lafleur et al., 2005a). This transition is generally rather abrupt, and is often associated with the acrotelm/catotelm transition in a peat profile (e.g., Rydin and Jeglum, 2006;Clymo, 1992). In HPM, the degree of decomposition is quantitatively tracked for each annual litter/peat cohort as the fraction of initial mass remaining (M/M o = i m i ÷ i m i,o ). Each PFT's factional contribution to a litter/peat cohort's mass can be quantified as m i /( i m i ). The degree of decomposition factor, M/M o or µ, is at the center of much of the dynamics and feedbacks in HPM. It influences decomposition rates (Eq. 7). Peat bulk density for litter/peat cohort j , ρ j , is modeled as a strongly non-linear function of µ j (Fig. 3) where ρ min is the minimum bulk density, ρ is the difference between minimum and maximum bulk densities, erfc is the complementary error function, c 3 and c 4 are parameters that control the rate of curvature of the erfc function and the value of µ j at which ρ j has increased half-way from minimum to maximum. Litter/peat bulk density affects (1) vertical peat accumulation rate, and thus water table depth (see Sect. 3.2) and vegetation productivity (Eq. 6); (2) peat water content in the unsaturated zone above the water table (see Sect. 3.2) and thus water table depth and peat decomposition rates (Eq. 7); and (3) peat hydraulic conductivity and thus run-off rates, peat water balance, and water table depth (see Sect. 3.2).

Water balance equations
Annual precipitation (P , m y −1 ) is a prescribed input variable. Evapotranspiration, ET, is calculated as a function of water table depth (Lafleur et al., 2005b;Ivanov, 1981), as where ET o is the maximum annual ET (m y −1 ), z 1 is the annual water table depth at which annual ET begins to decline, and z 2 is the annual water table depth at which annual ET reaches it minimum, and c 6 is the parameter needed to make a continuous ET function between z 1 and z 2 . Annual runoff, R (m y −1 ), is modeled as a function of annual water table depth and peat relative transmissivity, T , as where the first equation generates rapid drainage for flooded peatlands (i.e., z WT < 0), and R 1 is a base runoff rate proportional to peat depth, given by where R o sets the site-specific annual runoff (m y −1 ) for mean annual water table depth and peat depth of zero, c 8 is a parameter describing the rate of increase in annual runoff with increase in total peat depth, h PD . Relative transmissivity declines as the water table drops below the peat surface because the deeper, more decomposed peat has lower hydraulic conductivity. T is modeled as the ratio of the integrated hydraulic conductivity from the water table to the bottom of the peat to the integrated hydraulic conductivity for the entire profile where λ j is the thickness of litter/peat cohort j , T 0 is a parameter defining the minimum value of T , andK j the hydraulic conductivity of litter/peat cohort j , a function of bulk density (Walmsey 1977) is given by, The degree of saturation of litter/peat cohorts above the water table is determined by two factors -the cohort bulk density (ρ j ), which relates to pore size distribution and how readily peat dewaters under matric tension (e.g., Lafleur et al., 2005a), and the distance the cohort is above the water table (z WTz j , for z j < z WT ) and W min is the minimum peat water content (see Table 2). The absolute water content in the peat (for conservation of water) is the product of the degree of saturation, W , and cohort j peat porosity, j , where ρ om is the particle bulk density of organic matter. At the end of each time step after the water balance and carbon balance calculations have determined total peat height, peat cohort bulk densities and thickness, and total peat water content, the water table depth is determined by iterative search by annual cohort depth. The search identifies the water table depth for which the sum of water saturating all peat below that depth -i.e., ( j ·λ j ) for all cohorts below the water table -and the water occupying pore space in the unsaturated zone water content above that depth -i.e., (W j · j ·λ j ) for all cohorts above the water table -equals the total peat water content.

Model initialization and precipitation driver
HPM does not simulate peatland initialization explicitly. The model is initialized with a single fresh litter cohort (h PD ∼ 0.01 m) and z WT is held fixed at 0.07 m. These values of z WT and h PD determine vegetation composition and productivity (Eq. 6), and decomposition is calculated (Eq. 7), but the model does not calculate water table dynamics until h PD exceeds an initialization threshold (0.35 m) when the model initiates dynamic water balance calculations and z WT becomes a free variable.
A multi-millennial HPM simulation requires an annual reconstruction of the climate/weather. Such a record does not exist, but a regional Holocene climate reconstruction at 250year intervals has been developed from sediment, pollen and 8 S. Frolking et al.: A new model of Holocene peatland net primary production macrofossil records from the St. Lawrence lowlands and adjacent mountain areas (Muller et al., 2003b). Both the mean and confidence intervals of the Muller et al. (2003b) precipitation reconstruction were interpolated with a piecewise cubic Hermite interpolation (using Matlab function "interp1") to annual values, and truncated to generate an 8500-year Mer Bleue precipitation history, P M03 (t), and confidence interval, σ M03 (t). To add random variability to this record, we applied an autoregressive (AR[1]) model to generate random noise with persistence where ε is unit normal random number (and r(1) = ε), and φ (<1 for stationarity) determines the persistence. A random time series, r(1). . . r(8500), was normalized to a maximum absolute value of 1, r * (t), and then a precipitation time series to drive the model, P (t), was generated as where α is the amplitude of the maximum variance relative to σ M03 (see Fig. 4a; φ = 0.99, α = 2.5). In 1000 random P (t) series, ∼80% of the P (t) values fell within P M03 (t) ± σ M03 (t). We also conducted a simulation with constant precipitation over 8500 years at 0.94 m y −1 , the mean of P (t). Mean annual temperature in the region has been relatively stable for the past 8500 years (Muller et al., 2003b), so we have not included any temperature dynamics in the simulation.

Model evaluation data
To evaluate the HPM we sampled the variance in the longterm rates of carbon accumulation in the Mer Bleue peatland: a 28 km 2 cool temperate, raised peatland located ∼10 km east of Ottawa, Ontario (45.40 • N lat., 75.50 • W long., 69 m a.m.s.l.). Mer Bleue is a well-studied northern ombrotrophic bog (e.g., Roulet et al., 2007;Moore et al., 2002Moore et al., , 2007Lafleur et al., 2003Lafleur et al., , 2005aBubier et al., 2003Bubier et al., , 2006. The peatland formed over the past c. 8400 years, beginning as a fen and transitioning to a bog around 7100-6800 calendar years BP . Field measurements of contemporary carbon and water dynamics at the site include ∼10 years of eddy covariance net ecosystem exchange and surface energy balance measurements (e.g., Roulet et al., 2007;Lafleur et al., 2005a, b); chamber net ecosystem exchange of CO 2 (NEE) measurements in representative plant communities (e.g. Bubier et al., 2003); vegetation composition and biomass measurements (e.g., Moore et al., 2002;Bubier et al., 2006;Murphy et al., 2009); litter decomposition rate studies (e.g., Moore et al., 2007); and micrometeorological, water table, and run-off measurements. Mer Bleue is roughly oval shaped with an east-west orientation with two longitudinal lobes of fluvial sand/gravel material dissecting the western end of the bog, creating three separate arms (Fraser et al., 2001). We took a peat core a b d e c Fig. 4. (a) Four manifestations of annual precipitation (m y −1 ) with amplitude of random variability α = 2.5 and persistence φ = 0.99 (see . Solid black circles are regional mean annual precipitation at 250-year intervals, bars are regional variance, as determined by Muller et al. (2003b); solid black line represents constant mean annual precipitation. (b) HPM simulated annual water table depth (m below peat surface) for 4 random precipitation drivers and constant precipitation (solid black line). (c) Simulated annual total NPP for 4 random precipitation drivers and constant precipitation (solid black line). (d) Simulated annual total decomposition for 4 random precipitation drivers and constant precipitation (solid black line). (e) Simulated annual peat height for 4 random precipitation drivers and constant precipitation (solid black line). NOTE: all panels start with simulation year 50, avoiding the peatland initialization period (see text).  (Coûteaux, 1962), while the lower 4.53 m was sampled with a 1 m long, 0.075 m inner diameter Russian corer (Jowsey, 1966). C accumulation rates were established from an agedepth relationship combined with fine resolution measurements of density and C concentration in the cores.
At five depths the age was determined by radiocarbon dates (Table 3). Conventional radiocarbon dates are calibrated with CALIB 5.0 program (Stuiver and Reimer, 1993). The age for eight additional depths was determined from palynostratigraphical correlation with other radiocarbon dated pollen diagrams in the Montreal-Ottawa area (Muller et al., 2003a). In this we used BP to designate calendar date before present: present being 1950. Samples (1 cm 3 ) of fresh peat were dried at 105 • C, weighed, and ignited at 600 • C, to determine density and organic matter content, respectively (Dean, 1974;Beaudoin, 2003). Turunen et al. (2004) analyzed the top 0.60 m of peat for 48 cores from 24 different bogs in eastern Canada and obtained a mean of 0.46 g C g −1 , and this was used in our calculations. Measurements were made every 0.02 m from 0 to 2.00 m and every 0.04 m below that.
Long-term net rates of carbon accumulation were calculated using a smoothed, age-depth model including all age determinations for the core, the estimated C concentration for each sample, and the interpolated deposition time of each depth increment (e.g. either 0.02 or 0.04 m). The stratigraphy of the cores was described every 10 cm following the Troels-Smith (1955) nomenclature.

Sensitivity analysis
A series of simulations tested HPM sensitivity to model parameters; all were driven by stochastic precipitation parameterized as for the base runs, including several simulations that tested HPM sensitivity to the parameters controlling precipitation variability (φ and α in Eqs. 19, 20). Output variables examined include final peat height and carbon content, mean z WT during the final 40 years of the simulation, total NPP during the 8500-year simulation, and the percent of that total NPP remaining as peat at the end of the simulation.

Results
HPM simulation output includes both the time series of annual net carbon and water fluxes and peat accumulation, as well as a developing peat profile. At the end of a multi-millennia simulation, the simulated peat profile can be "cored" to get simulated contemporary profiles of peat properties (depth, age, bulk density, moss and vascular plant fractions of remaining peat) that can be compared to field cores.

Model evaluation against site data at Mer Bleue Bog
After the initialization period (∼25 years in this scenario), the annual water balance is calculated, the resulting water table depth is determined (conserving total water), an annual carbon balance is calculated (productivity minus decomposition), and a final water table depth is determined (again conserving water). The simulated water table rapidly drops from the initiation value of 0.07 m to about 0.25 m once the water table is allowed to be dynamic. For most of the simulation the water table fluctuates between about 0.2 m below the peat surface during the wet periods and 0.35 m below the peat surface, dropping below 0.4 m during dry periods (Fig. 4b). The magnitude of interannual variability in z WT was higher during periods of low precipitation and deeper water tables (Fig. 4b). Total NPP drops from initial rates of about 1.4 kg C m −2 y −1 to about 0.7 kg m −2 y −1 after about 1000 years, due to the transition from minerotrophic to ombrotrophic PFT dominance, and then declines more slowly to about 0.43 ± 0.01 kg C m −2 y −1 during the final millennium of the simulation (Fig. 4c). Decomposition also drops from initial rates of about 1.2 kg C m −2 y −1 to about 0.7 kg C m −2 y −1 after about 1000 years, and then declines more slowly to about 0.41±0.01 kg C m −2 y −1 during the final millennium of the simulation (Fig. 4d). Decomposition is less sensitive to wet and dry periods than NPP, but increases during periods of declining precipitation and increasing water table depth. Simulated total peat height rises to about 1.5 m in the first 1000 years, then more slowly to about 4 m after 5000 years, then irregularly rises and falls for about 1500 years before rising again more slowly to about 5.4 m at the end of the 8500 years simulation (Fig. 4e).
The observed age-depth profile in the MB930 core indicates relatively rapid peat accumulation prior to ca. 5500 BP, an interval of slow accumulation ca. 5500 BP-ca. 2500 BP, and then more rapid accumulation again during the past ∼2500 years (Fig. 5a). Total simulated NPP over the 8500year simulation was about 4800 kg C m −2 , with 6.5% of that remaining in the final peat core (310 kg C m −2 ; Table 4). Mean annual NPP was about 0.56 kg C m −2 y −1 , and mean long-term peat accumulation was about 0.037 kg C m −2 y −1 . Total minerotrophic PFT NPP was about 35% of the total, but accounted for only about 20% of the peat; total vascular PFT NPP was also about 65% of the total, but accounted a b c annual total NECB (kg C m -2 y -1 ) for only about 35% of the peat (Table 4). Only the brown and Sphagnum moss PFTs had more than 10% of their total NPP remaining in the peat at the end of the simulation ( Table 4). The mean accumulation rate over 8500 years was about 0.6 mm y −1 and 35 g C m −2 y −1 . The simulated peat is younger at all depths except ∼2 m than observed MB930 core dates (Fig. 5a). The discrepancies between the simulated and observed profile are largest in the most recent 1000 years (∼1 m of peat in the simulation, but ∼0.5 m in the peat core) and in peat 3800-5600 years old (∼1.1 m thick in the simulation and ∼0.5 m thick in the peat core) (Fig. 5a).
The net rate of C accumulation (NRCA), estimated from the core as the carbon content per cm depth divided by the age accumulation per cm depth, ranges from ∼0.01 kg C m −2 y −1 to about 0.07 kg C m −2 y −1 . The ob-served NRCA had high and highly variable rates before 7000 BP, a small increase around 6000 BP, then a low rate 5500-3000 BP, rising to a moderate and moderately variable rate between 3000 BP and the top of the peat core (Fig. 5b). In the simulation, each annual cohort's mass at the end of the simulation is equivalent to reported NRCA values estimated from dated cores, and ranged from ∼0.015 g C m −2 y −1 to ∼0.08 g C m −2 y −1 (Fig. 5b). The simulated pattern is generally similar, though with several differences. After the rise in apparent accumulation rate around 6000 BP, the simulated NRCA drops more steadily to a low value around 3000 BP, and then increases less abruptly. As in the age depth profile (Fig. 5a), the simulated carbon accumulation rate for the shallow peat is significantly larger than in the observed core. Although the final core age-depth profile shows monotonically increasing depth with age (Fig. 5a), the simulated peat height does not increase monotonically (Fig. 4e). The net ecosystem annual carbon balance (NECB) of the peatland (total NPP minus total decomposition in HPM) varies between a net uptake of about 0.2 kg C m 2 y −1 and a net loss of −0.15 kg C m 2 y −1 (Fig. 5c). The variability is driven by precipitation variability and resulting variability in water table depth. In general, and particularly during the simulation years 5000-7000, the variability in NECB results mostly from variability in productivity (Fig. 4c), not variability in decomposition (Fig. 4d).
Since HPM calculates NPP and litter/peat decomposition by PFT, the simulated core records the composition of the remaining peat by PFT and this can be compared with the results of observed PFT remains deduced from macrofossil analysis. The macrofossil analysis indicates the deeper peat is predominantly herbaceous, overlain by a meter or so with a significant woody peat fraction, and bryophyte peat dominating the upper 2.5 m. Throughout the core a significant fraction of the peat consists of Detritus granosus, which constitute the partially decomposed peat matrix. The state of humification of the entire material varies throughout the core (Fig. 6c). In all four base-run simulations the older peat (or deeper core) is dominated by minerotrophic PFT peat (roughly 7000-8500 BP) and the younger peat is dominated by ombrotrophic PFT peat (Fig. 6a, b). For most of the core, lawn and hummock Sphagnum peat comprises more than half the total; the herbaceous fraction was greatest deep in the core, and the woody fraction averaged around 15% (Fig. 6a, d).
The simulation with constant annual precipitation for 8500 years generated a smoothly monotonic water table depth dropping from 0.20 to 0.35 m below the surface, a smoothly monotonic and declining NPP, decomposition, and NECB, a smoothly monotonic peat height curve rising to a final height of 4.4 m (and a final carbon content of 250 kg C m −2 ), a smooth age-depth profile concave upward from present back to about 6000 BP (bog phase) then concave downward back to 8500 BP, and a smooth NRCA curve with a broad minimum from about 4000-1000 BP (see . The constant precipitation simulation had more rapid accumulation than the base runs during the first 2000-3000 years of the simulation (Fig. 4e), despite lower NPP (Fig. 4d); during this period the constant precipitation rate was greater than the Muller et al. (2003b) reconstruction (Fig. 4a), resulting in a generally shallower water table (Fig. 4b) and lower decomposition (Fig. 4d).

Peat bulk density parameters
Reducing ρ min from 50 to 35 kg m −3 ( ρ remained 70 kg m −3 , so maximum bulk density was also reduced by 15 kg m −3 ) had little impact on the total peat carbon accumulation, but increased the final peat height by about 1 m (∼20%) and resulted in a ∼0.08 m deeper water table (SR#1 in Table 5). Shallow, relatively undecomposed peat has lower bulk density, so in this simulation each litter cohort was thicker (same mass, lower ρ), but decomposed at approximately the same rate, so took about the same time to "collapse" to higher bulk density. In HPM, the water table generally is located near this bulk density transition (unless it is an extreme precipitation year) so the mean water table was deeper with lower ρ min . In addition, runoff increases with peat height (Eq. 13), so in HPM deeper peat has higher runoff and a deeper water table, if hydraulic transmissivity can accommodate it. Despite a deeper water table, the fraction of total NPP remaining as peat increased to 7.8% due to a higher fraction of NPP by ombrotrophic PFTs that prefer a deeper water table and have lower initial decomposition rates (Table 1). Reducing ρ from 70 to 50 kg m −3 increased the final water table depth by about 0.04 m, reduced both total NPP and the percent remaining as peat, so total peat carbon accumulation decreased by ∼20%, though total peat depth by only ∼5% due to decreased overall bulk density (SR#2). The reduction in peat depth was less than the reduction in peat mass because the peat column overall bulk density was decreased. Increasing ρ from 70 to 90 kg m −3 increased total NPP and the percent remaining as peat and thus total peat mass by ∼13%, but total peat height was slightly lower due to increased overall bulk density (Table 5). Parameterizing the peat collapse to start earlier and/or proceed more slowly (SR#4-6, Table 5, Fig. 3) reduced peat accumulation and the fraction of NPP remaining as peat and caused a shallower final water table. A more abrupt transition (SR#7) had little effect.

Vegetation productivity and decomposition rates
Reducing NPP for all PFTs by 25% reduced total peat mass and peat height by ∼25%, and resulted in a ∼0.1 m shallower water table (SR#8 in Table 5). A shallower water table would result from both lower runoff due to less peat height and a shallower transition to high bulk density, low transmissivity peat due to thinner annual cohorts. Note that total NPP was only reduced by 12.5%, in part because the peat height took ∼1000 years longer to reach the transition to ombrotrophy and lower relative NPP (though difference in water table would also have an effect). However, the fraction remaining was also reduced from 6.5% to 5.4%, because a greater fraction of total NPP came from relatively more decomposable minerotrophic PFTs. Increasing NPP for all PFTs by 25% had the opposite effect (SR#9). Similarly, increasing the initial decomposition rate for all PFTs by 25% decreased total peat mass and peat height (by ∼10%), and resulted in a 0.07 m shallower water table (SR#11). Again, a shallower water table would result from both lower runoff potential due to less peat height and a shallower transition to high bulk density, low transmissivity peat due to thinner annual cohorts. Total NPP increased by 8%, due to a slower transition to ombrotrophy, but only 5.4% remained as peat, despite  a shallower water table, due to higher initial decomposition rates. Decreasing initial decomposition rates had an opposite and slightly enhanced effect (SR#10). The anoxia scale length parameter (c 2 , Eq. 9, Table 2) controls the decline in overall decomposition rate of submerged peat as a function of depth below the water table. Increasing the scale length value to 0.5 m resulted in a ∼25% decrease in peat depth and mass after 8500 years, and 0.04 m reduction in final water table depth; decreasing the value to 0.2 m resulted in a ∼10% increase in peat depth and mass after 8500 years, and no significant change in final water table depth (SR#12-13 in Table 5).
In a simulation where all vascular NPP was input at the peat surface (i.e., same total NPP, but no roots), the final peat accumulation was reduced by only a few percent (both total depth and total mass; SR#14 in Table 5), which may not be significant, given the stochastic factor in the precipitation. The water table depth in this simulation was ∼0.05 m shallower than in the base runs; this was likely the result of the peat at ∼0.2-0.25 m being more decomposed (higher bulk density) in this simulation (no fresh litter inputs below the surface), and thus lower transmissivity (Eq. 14) reduced runoff. This would reduce overall decomposition slightly, apparently nearly offsetting the increased overall decomposition of 100% surface litter, and minimizing the impact on total peat accumulation. Reducing sedge PFT belowground NPP fraction from 0.8 to 0.5 had little impact (SR#15 in Table 5).

Water balance parameters
Varying the amplitude of the precipitation variability (α) from 1.0 to 3.0 had little impact on the results (SR#16-19; Fig. 7). Reducing the persistence parameter for the random precipitation (φ; base value = 0.99) had little impact on total peat height or depth or final water table (SR#20-23, Fig. 8), but total NPP and percent remaining as peat were different (with offsetting changes in terms of final peat carbon content). Increasing the persistence parameter for the random precipitation caused substantial variability in final peat height, peat depth, and final water table (SR#24-27, Fig. 9b) as well as in the dynamics of peat accumulation and its representation in the final peat core age-depth profile and carbon density per unit time (Figs. 8c-e, 9c-e).
Since precipitation is based on a paleo-reconstruction (Muller et al., 2003b), evapotranspiration is parameterized from contemporary eddy covariance measurements (Lafleur et al., 2005b), and run-on is assumed to be negligible, the "free" variable in the simulated annual water balance is annual run-off. Three parameters control the runoff -R 0 determines the base annual run-off from P-ET 0 for a shallow, saturated peat; c 8 controls the linear rate of increase in runoff with peat height, and T 0 sets a minimum value for deep water table reduction in annual transmissivity relative to the rate for saturated peat (Table 2 and Eqs. 13,14). All results were sensitive to all three parameters -in particular, varying these parameters caused the greatest variability in the percent of total NPP remaining as peat (SR#28-39 in Table 5). All run-off parameters affected the water table depth, with a clear pattern of shallower (deeper) water tables correlating with a greater (smaller) percentage of total NPP remaining a peat than in the base runs. For example, reducing c 8 from 0.2 to 0.3 caused a 0.06 m decrease in z WT , and an increase in final peat height and carbon content by almost 35%, despite a 7% reduction in total simulation NPP, as ∼9.3% of that NPP did not decompose, compared to 6.5% in the base runs.

Discussion
HPM simulates the long-term carbon and water dynamics of a peatland, calculating annual carbon and water balances for each year of a multi-millennium simulation. Most of these results are not possible to test directly; there are multiyear or annual carbon balance budgets for a few peatlands (e.g., Roulet et al., 2007;Nilsson et al., 2008;Koehler et al., 2010), but these are on peatlands that are at least several thousand years old, and they may not be representative of the carbon balance of those peatlands throughout their development. Contemporary C balance records of 5-10 years put clear constraints on the current system state and dynamics. However, they provide no information on how these systems behaved in the past (vegetation dynamics, carbon accumulation, disturbance/response) and little information on how these systems will behave under climate/weather conditions outside the range of variability over the past 5-10 years. Contemporary carbon and vegetation dynamics measurements from a set of peatlands representing a chronosequence of stages of development (i.e. space for time substitution) within the same climatic region (e.g., Leppälä et al., 2008;Merilä et al., 2006) can provide a series of snapshots of peatland behavior through time.
Most information on past dynamics comes from peat core analyses; data include peat age, carbon content, bulk density, macrofossil identification and pollen content, and testate amoebae community composition. This information is then used to generate estimates of the dynamics of peat vegetation composition (e.g., Mauquoy et al., 2008), past climatic conditions (e.g., Mauquoy et al., 2008) and water table dynamics (e.g., Charman, 2007), and long-term C accumulation rates (e.g., Turunen, 2003) -all of which are usually explained inductively. HPM generates a number of variables that are comparable to these contemporary observables -final state (i.e., contemporary) peat depth and total mass or carbon content, and vertical profiles (at a resolution of 1 year or <0.01 m) of peat age, bulk density, degree of decomposition, and PFT composition. HPM results are generally consistent with long-term peat core records from Mer Bleue, i.e., about 5 m and 300 kg m −2 peat carbon, apparently accumulating at variable rates over the past 8500 years, with deeper (older) peat dominated by minerotrophic PFTs and shallower peat dominated by ombrotrophic PFTs, though the composition of the simulated core seems to be less variable than observed in the field (Figs. 5 and 6). Bryophyte remains account for ∼65% of the total peat (Table 4); Turetsky (2003) reported that bryophytes are estimated to comprise about half of the total peat, based on more than 600 archived peat cores from continental Canada cataloged by Zoltai et al. (2000). This result is also consistent with the observations of the importance of Sphagnum in determining shifts in periods of peatland growth (Belyea and Malmer, 2004).

Vegetation
Due to large variability in northern peatlands, very few HPM model parameters (Tables 1 and 2) are well constrained by field data. Primary productivity is highly variable between sites and in different years for Sphagnum mosses, true mosses, forbs, shrubs, and trees (e.g., Wieder, 2006, and references therein) and for total site NPP (e.g., Rydin and Jeglum, 2006;Moore et al., 2002, and references therein). Relative NPP for the PFTs in HPM, and their sensitivities to water table depth and peat height are meant to represent general patterns common to northern peatlands (e.g., Walker, 1970;Reader and Stewart, 1972;Bernard, 1974;Forrest and Smith, 1975;Tallis, 1983;Backeus, 1990;Korhola, 1992;Rydin, 1993;Klinger and Short, 1996;Hughes and Barber, 2004;Leppälä, et al., 2008;Murphy et al., 2009;Murphy, 2009 and references therein). It is important to note that model NPP for a given water table and peat depth combines both plant abundances under those conditions as well as PFT inherent productivity per unit leaf biomass, so the defined optima and tolerance along water table and peat depth gradients represent PFT ecological niches instead of wider physiological niches. For example, while ombrotrophic sedges could be relatively productive with shallow peat, they are likely to be outcompeted for space by minerotrophic species and thus are parameterized to have low NPP.
This representation of relative NPP results in a simple dynamic vegetation model, but does not simulate PFT establishment and associated lags; the assumption is that at least one representative species of each PFTs is effectively present, even if with essentially zero productivity, and becomes productive when conditions (water table depth and peat depth) are suitable; common peatland plants are able to disperse well via air and water (Salonen, 1987;Talbot, 2009). The transition from minerotrophic PFTs to predominantly ombrotrophic PFTs (fen-bog transition) is parameterized to occur as the peat depth increases from roughly 1 to 2 m; this transition depth is likely dependent in a site-specific way on the surface-and groundwater hydrology and miner-alogy of the watershed encompassing the peatland (Fraser et al., 2001;Comas et al., 2004;Glaser et al., 1990). Macrofossil, pollen, and radiocarbon analyses indicates that Mer Bleue transitioned from a fen to a bog circa 7100-6800 BP; in the base-run simulations the minerotrophic PFT fraction of contemporary peat drops from ∼40% to ∼25% over the age range of 7500-6500 BP (Fig. 6). HPM results do not seem to be very sensitive to these above-and below-ground allocation fractions (Table 5).
HPM uses 12 PFTs to represent the range of vegetation found on many peatlands. This collection of PFTs was determined by first thinking of key functional differences: aboveand below-ground biomass/NPP; vascular or non-vascular; ombrotrophic or minerotrophic; woody or herbaceous; tissue decomposability. Ombrotrophic mosses were then partitioned into a range of "moisture" classes (hollow, lawn, hummock Sphagna and feathermoss), as these are often used in paleo-analysis as a wetness indicator. Differences in productivity (e.g., Gunnarsson, 2005;Pedersen, 1975) and decomposability (e.g., Hogg, 1993;Belyea, 1996) will influence rates of carbon accumulation as peat. One missing PFT is trees (unless trees can be considered to be just large shrubs); trees are likely to be important for simulating many current peatlands (particularly in the tropics), and also to capture the impacts of climate drying on peatlands carbon and water dynamics. We believe that HPM is a first attempt to quantitatively link long-term peat carbon accumulation and hydrology to vegetation composition and dynamics in a process-based model. With 12 PFTs, HPM could be used to look at the response of individual plant functional types to past and future climatic changes. While 12 (or more) wetland/peatland PFTs may seem too few from an peatland ecology point of view, it is likely to be more vegetation classes than a global earth system model will be willing to incorporate; for example, Wania et al. (2009) used only two wetland PFTs in their LPJ-Why model -flood-tolerant C3 graminoids, and Sphagnum moss. A series of simulations with HPM would provide a mechanism for evaluating the importance of different PFTs in long-term peatland dynamics and carbon cycling.

Decomposition
Decomposition rates are generally low in boreal peatlands, and plant tissue quality contributes to this through low nutrient content and high refractory content (Moore and Basiliko, 2006;Laiho et al., 2006). In situ litter-bag decomposition studies show general patterns in initial decomposition rates among plant functional types (forbs ≥ sedges > shrubs > Sphagna > wood) but within-PFT variability is high and PFT ranges generally overlap (e.g., Moore et al., 2007). Litter bag decomposition rates are determined by fitting mass loss data to an exponential function, m(t) = m 0 e −kt , under the assumption that k is a constant (dm/dt=-km); since HPM represents litter/peat decomposition as an initial rate k 0 , that declines linearly with mass loss (i.e., dm/dt = −(m/m 0 ) k 0 m), the litter bag k-values need to be modified to k 0 values. We have modified these so that, after 5 years of decomposition they have approximately the same mass loss (see Table 1).
The anoxia scale length parameter (c 2 , Eq. 9, Table 2) controls the decline in overall decomposition rate of submerged peat as a function of depth below the water table. In situ biogeochemical studies indicate that decomposition rates in deep peat are limited by lack of solute transport (e.g., Beer and Blodau, 2007;Beer et al., 2008). In HPM this is represented as a rate multiplier that declines exponentially with distance below the simulated water table. The base run value for this scale length is 0.3 m, which is roughly the range of interannual variability (∼0.35 to 0.70 m) in maximum growing season water table depth observed at Mer Bleue . The minimum decay rate multiplier is set to 0.001 y −1 -this does not include reduction in decomposition rate from the initial value due to mass loss (see Eq. 7; at depth total cohort m/m 0 < 0.1). Simulation results were not sensitive to this minimum value. Increasing this scale length means that "full anoxia" is pushed deeper into the submerged peat; this could be caused by larger seasonal to interannual variability in the water table, or by non-negligible water (and solute) flow through the submerged peat. Increasing the scale length value to 0.5 m resulted in a ∼25% decrease in peat depth and mass after 8500 years, and 0.04 m reduction in final water table depth; decreasing the value to 0.2 m resulted in a ∼10% increase in peat depth and mass after 8500 years, and no significant change in final water table depth (Table 5).

Hydrology
Variability in peat net carbon uptake and in the contemporary carbon density per unit time down the profile is driven in HPM by variability in precipitation that drives a variability in water table depth, providing the entry in to the interconnected carbon and water cycles represented in Fig. 1. Constant precipitation generates a smoothly varying peat core (Figs. 5b,6b). Increasing the amplitude of precipitation variability causes high frequency variability in carbon and water variables (e.g., Fig. 7), but only when these deviations from the mean are persistent over several decades to centuries are there pronounced deviations in peat height or the carbon density or age depth profiles (Figs. 8c-e, 9c-e). Since precipitation and snowmelt inputs and runoff are highly episodic and evapotranspiration is strongly seasonal, peatland water tables are seasonally dynamic, and data are often missing or difficult to interpret in winter due to freezing conditions (e.g., Lafleur et al., 2005b. HPM has an annual time step with no seasonality in water inputs or outputs, and computes a mean annual water table, and it is not clear what exactly this should be compared to. Does it represent an observed actual mean annual water table, or a mean unfrozen season water table, or neither? Should its absolute magnitude be compared to field observations, or are relative interannual variability and long-term trends more important? Charman (2007) showed that proxies for peatland water table reconstructions are most likely related to summer water deficits, which themselves are most related to summer precipitation and, to a lesser degree, summer temperature (evaporative demand). Charman (2007) suggests that lowfrequency variability in annual moisture deficits (i.e., decadal or longer), leading to persistent wet or dry conditions, can cause changes in plant community composition and decay rates, which feedback, through changes in peat hydraulic properties, on peat surface wetness.

Net carbon balance
Total decomposition must be less than total NPP over time in order for peat to accumulate. In the 5 base-run simulations total decomposition averaged 93.5% of total NPP (range: 93.4%-93.7%; Table 5); total NPP had a standard deviation of 57 kg C m −2 (1.2% of the mean) and total decomposition 58 kg C m −2 (1.3% of the mean), but total peat mass had a standard deviation of only 4 kg C m −2 (also 1.3% of its mean), implying that NPP and decomposition co-vary. Across all simulations in Table 5, total decomposition was highly correlated with total NPP (decomposition = 1.05 · NPP -565; R 2 = 0.99), and total decomposition ranged from 90.5% to 96.1% of total NPP. With total NPP ∼5000 kg C m −2 , a change in fractional decomposition from 94% to 93% implies 50 kg C m −2 less peat in the final profile, about 17% of the total, so net peat carbon accumulation is the relatively small difference of two much larger gross carbon fluxes (NPP and decomposition). On a much shorter timescale, Bubier et al. (1998) found a strong correlation between chamber measurements of gross photosynthesis and ecosystem respiration (as CO 2 ) across a range of boreal Canadian peatlands from rich fen to bog.
Peatland NRCA is measured on contemporary cores, and must always be positive, as one cannot sample lost carbon in a core but only infer loss, either from the presence of charcoal or from long time intervals separated by little vertical distance. Peatland NECB is annual change in storage, and this can be positive or negative, but only can be quantified by direct measurement during the year in question. Simulated NECB shifted from a positive to a negative carbon balance throughout the simulation; as a result, the correlation between simulated NECB and simulated NRCA was only moderate (r 2 ∼0.4-0.5). NECB variability results mostly from variability in productivity (Fig. 4c), not variability in decomposition (Fig. 4d), consistent with contemporary field studies (e.g., Riutta et al., 2007;Laine et al., 2009a). Interannual variability in NECB (Fig. 5c) is much less evident in the final peat core age-depth profile or field NRCA (Fig. 5a, b) because millennia of decomposition reduce the magnitude of variability and field sampling inevitably aggregates over a number of years, acting as a low-pass filter. Note that HPM simulates total decomposition (mass loss), but does not disaggregate this into CO 2 , CH 4 , and DOC Nilsson et al., 2008;Koehler et al., 2010). Based on the HPM simulations, we hypothesize that peatlands can have a weak negative C balance for decades or even centuries (e.g., Figs. 4e, 5c, 7c, 8c, 9c) that cannot be unambiguously determined from dating a peat core profile. These intervals do not need to have low NPP, and so the peatland may not be severely stressed.

Outstanding issues and future directions
This first version of HPM does not yet include several factors that will be important to develop for more widespread application. Mean annual temperature exerts controls on both NPP and decomposition; we anticipate that, as with HPM sensitivity to water table depth, long-term variability and trends (e.g., Holocene Thermal Maximum, Medieval Warm Period, Little Ice Age) will play a more important role than short-term interannual variability. Temperature effects become particularly dramatic when they generate or thaw permafrost (e.g., Robinson and Moore, 2000;Turetsky et al., 2002a). Trends in seasonality in temperature (e.g., related to changes in summer insolation; Jones and Yu, 2010) and growing season length may also be important. In their regional paleoclimate analysis, Muller et al. (2003b) detected a long-term trend towards more frequent cooler summers and more frequent warmer winters (i.e., decreasing seasonal contrast). Seasonal hydrological dynamics (e.g., summer droughts in years with normal annual precipitation) could have a large impact (e.g., Charman, 2007), but will be difficult to reconstruct for past millennia. Peatland disturbances, particularly fire, but more recently human impacts, are not simulated. Mer Bleue core MB930 had evidence of charcoal at about 1.8 m and 4.6 m. If these fires severely burned the peat (e.g., Zoltai et al., 1998;Robinson and Moore, 2000;Turetsky et al., 2002b), they would have a significant impact on the overall peat accumulation and age-depth profile.
Several other factors are included in a simple way that may preclude additional feedback dynamics in peatland development. For example, nutrient and mineral input and mineralization dynamics can influence system behavior (e.g., Pastor et al., 2002), but in HPM are represented only through a peat depth proxy. Vegetation dynamics are represented through the smooth functions of Eq. 6 (see Fig. 2), with no stochasticity, stress threshold response, or dispersal lags, which are likely to make any particular site more variable than the model output.

Conclusions
In essence, HPM is a set of dynamically linked multiple hypotheses on the complex nature of a strong coupling between peatland annual carbon and water balances (Fig. 1). The linkages in HPM are unambiguous and therefore in princi-ple could be deductively tested (i.e., subject to falsification) if appropriate observations exist, but the conceptual framework, links and general structure is based on the inductive inferences of peatland development that shape how we think peatlands operate as a system (e.g. Charman, 2002;Rydin and Jeglum, 2006). The links between precipitation, water table, vegetation composition and productivity, decomposition, and peat hydraulic properties generate a relatively stable system overall. Variable precipitation leads to variability in simulated water table depth, which leads to variability in the carbon balance. This is most apparent in the final core mass as a function of age (Figs. 5b, 7d, 8d, 9d), which integrates the combined effects of relative productivity and total decomposition, both dependent on water table dynamics over years (productivity) to decades or centuries (decomposition). Variability in net carbon accumulation that would be measurable in a peat core (e.g., Fig. 5b), and that might have an impact on long-term carbon accumulation, results from persistent (tens to hundreds of years) precipitation deviations or trends (order of 10%) away from a long-term mean (e.g., Figs. 8 and 9). Higher frequency variability gets smoothed out in the long-term record by the "low-pass filter" of decomposition and sampling.