Abstract
Northern peatlands are a globally significant carbon (C) reservoir, yet also function as dynamic methane (CH4) sources to the atmosphere. The fate of peatland C stores and related climate system feedbacks remain uncertain under scenarios of a changing climate and enhanced anthropogenic pressure. Here, we present a synthesis of Holocene peatland C dynamics for the Hudson Bay Lowlands (HBL), Canada, in relation to the past atmospheric CH4 trends, glacial isostatic adjustment, and paleoclimate. We report that peatland age and trophic status, together with paleoclimate, contribute to explaining some of the temporal variation in C accumulation rates (CARs) in the HBL. Our results show that younger, minerotrophic peatlands accumulate C faster, and although detailed paleoclimate data are not available, the results suggest the possibility of higher CARs in association with warmer Holocene climates. Peat initiation rates and CARs were greatest during the mid-Holocene; however, our results reveal that two-thirds of the HBL C pool is stored in peat of late Holocene age, owing to long-term peatland expansion and development. Whereas the HBL has been a net C sink since mid-Holocene peat initiation, the HBL also appears to have been a modest C source, with 85% of the losses occurring during the late Holocene as a consequence of the gradual decay of previously accrued peat. Late Holocene peat decay, under wetter climatic conditions, and from a landscape occupied by an abundance of minerotrophic peatlands, indicates that the HBL may have been a natural terrestrial source of CH4 to the late Holocene atmosphere. While the peatlands of the HBL may continue to function as a globally significant C store, ongoing C losses from the HBL may have important implications for the global C budget and climate system.
Keywords
Introduction
Northern peatlands are important repositories of atmospheric carbon (C), storing at least 500 Pg C or about one-third of the global C pool (Yu, 2011), within approximately 3% of the terrestrial surface area (Clymo et al., 1998). Although global syntheses confirm the importance of peatlands as net C sinks in the global C cycle and as cooling agents in the global climate system (Frolking et al., 2011; MacDonald et al., 2006; Yu, 2011), peatlands also function as important C-sources, principally in the form of carbon dioxide (CO2), methane (CH4), and dissolved organic C (DOC). Anticipated high-latitude changes in temperature and the net moisture balance may enhance both peatland net primary productivity (NPP) and decomposition. Recent evidence using peat records covering the last millennium suggests that NPP may be more important than peat decay in determining the rate of C sequestration (Charman et al., 2013; Yu, 2012). Furthermore, peat initiation and long-term development appear to be sensitive to climate; however, the relationship between paleoclimate and temporal variability in the rate of peat C accumulation that is chronologically controlled using multiple vertical peat dates (referred to here as C accumulation rate (CAR)), remains poorly constrained (Charman et al., 2013). As a result, improving our understanding of mechanisms that could augment peat C decomposition and CO2/CH4 emissions, which together could result in a positive feedback to climatic warming, continues to be a high priority (Jones and Yu, 2010). Assessments of long-term peatland-climate-C dynamics provide an important perspective for informing our understanding of the future trajectory of peatland C accumulation, under conditions of a changing climate and enhanced anthropogenic pressure.
Long-term peatland development is influenced by the interaction of autogenic and allogenic processes and associated feedbacks. As a consequence of an imbalance between productivity and decay under generally waterlogged conditions, peat C accumulates over millennia. Examination of long-term peatland-climate-C cycle dynamics relies upon assessments of the synchronicity between regional to global syntheses of peat initiation and expansion, trophic status, and CAR in the context of inferred paleoclimate and long-term atmospheric trace gas variation. Results from these approaches provide evidence that the majority of northern peatland initiation occurred in the early Holocene, under conditions of maximum summer insolation and temperature seasonality (Yu et al., 2010), synchronous with a rapid rise in ice core inferred, atmospheric CH4 concentrations (MacDonald et al., 2006). Furthermore, the development of CH4 emitting pre-peatland marsh conditions coupled with rapid peatland initiation and expansion of minerotrophic systems (both sources of CH4) has been linked to the early Holocene atmospheric CH4 rise (Jones and Yu, 2010; Smith et al., 2004; Yu et al., 2013).
Late Holocene polar ice-core records reveal a second, more gradual rise of atmospheric CH4 concentrations, following a mid-Holocene reduction. However, more uncertainty exists regarding the contribution of northern peatlands to this rise in atmospheric CH4 during the mid- and late Holocene. Apparent declines in the rate of northern peatland expansion, coupled with potential ombrotrophication (fen to bog transition) of existing peatlands, may contribute to reductions in peatland-derived CH4 emissions (Yu, 2011). However, recent evidence from the extensive peatland region of northeastern Canada, document delayed ombrotrophication until the latter half of the late Holocene (Holmquist and MacDonald, 2014), and a late Holocene, possibly climate-driven, return of minerotrophic conditions (Van Bellen et al., 2013). Both scenarios would favor potential late Holocene, peatland-derived CH4 emissions. However, late Holocene permafrost establishment (Kuhry, 1998, 2008; Lamarre et al., 2012) may contribute to local reductions in potential CH4 emissions, which may be remobilized as established permafrost degrades under climate warming or disturbance scenarios (Kuhry et al., 2010).
Alternative hypotheses assert that possible terrestrial CH4-sources during the mid- to late Holocene CH4 rise may be related to early anthropogenic activity, such as rice cultivation and deforestation (Ruddiman, 2003), tropical wetland contributions (Brook et al., 2000; Yu et al., 2010), and/or lateral expansion of existing high-latitude peatlands through the mid- to late Holocene (Korhola et al., 2010). Until recently (Holmquist and MacDonald, 2014; Packalen et al., 2014), potential late Holocene CH4 emission estimates from northern peatlands have been limited by a lack of spatially explicit evidence of mid- to late Holocene peatland initiation and expansion dynamics from major peatland regions, such as the Hudson Bay Lowlands (HBL), Canada. In consideration of new evidence, and building upon the hypothesis of potential CH4 contributions related to peatland expansion (Korhola et al., 2010), together with the possibility of sustained minero- to oligotrophic conditions in important northern peatland regions, we suggest that further analyses may indicate some CH4 contributions from northern peatlands to the late Holocene CH4 rise (Behl, 2011; Packalen et al., 2014; Yu, 2012).
The timing of peat initiation in the HBL appears to be principally controlled by glacial isostatic adjustment (GIA) following retreat of the Laurentide ice sheet and land emergence from the post-glacial Tyrrell Sea (Glaser et al., 2004a, 2004b; Packalen et al., 2014). As well, evidence of peat initiation periodicity reveals that the most rapid and intense period of peat initiation in the HBL occurred during the mid-Holocene, in advance of the late Holocene CH4 rise (Packalen et al., 2014). Although peatland succession trajectories in the HBL can be diverse (Klinger et al., 1994; Sjörs, 1959), newly initiated peatlands begin as nutrient-rich, minerotrophic systems that are dominated by herbaceous plant communities and more easily decomposable plant litter. Following millennia of succession, peatlands may transition to more nutrient-poor, ombrotrophic systems dominated by mossy vegetation producing litter that is more difficult to decompose. The trajectory between these two endpoints contributes to variation in long-term peatland C cycling.
While evidence suggests that climate is an important control on northern peatland initiation and C accumulation dynamics, the link between regional climate in the HBL and variation in the rate of C accumulation is not well documented (Gorham et al., 2012). Delayed glacial retreat in northeastern Canada likely deferred the Holocene Thermal Maximum (HTM) until the mid-Holocene, as orbitally driven insolation and seasonal differences in temperature were declining (Renssen et al., 2009). Millennial-scale climate variation in the vicinity of the HBL is approximated by a synthesis of pollen-inferred temperature and precipitation records obtained from northern Quebec, Canada (Viau and Gajewski, 2009). Although northern Quebec and the HBL are not identical in biophysical or physiographic environments, both share similar latitudes and are influenced by the Hudson Bay. These pollen-based temperature reconstructions indicate little temperature change since initiation of HBL peatlands began ~8 cal. kyr BP. Reconstructed Holocene temperature fluctuations are in the order of ±1°C relative to modern climate (1960–1990 climate normal), with mostly cooler temperatures during the mid- and late Holocene, and evidence of similar summer temperatures and a trend toward warmer winters relative to today between 4.5 and 2 cal.kyr BP. Reconstructed precipitation changes show a shift at about 3 cal. kyr BP from somewhat reduced precipitation in the earlier part of the record to increased precipitation in the late Holocene (Viau and Gajewski, 2009).
Here, we present the variability in Holocene C balance in relation to the key potential drivers of that variability, including major temporal trends in paleoclimate, ice-core inferred trace gas variation, and frequency of regional peat initiation events in the HBL. Working under the assumption that both allogenic and autogenic processes account for variation in the rate of C sequestration in the HBL, we quantify CAR and estimate long-term decay in the HBL, at 0.5-kyr intervals to model changes in the net C balance (NCB) through time. Using our analysis of CAR in the HBL, together with a synthesis of all available peat initiation records for the HBL (Packalen et al., 2014), we model net C uptake (NCU) and release terms at 0.5-kyr intervals following the method of Yu (2011). We then compare the resulting patterns in C uptake and release in the HBL with mid- to late Holocene, pollen-inferred paleoclimate. We hypothesize that reduced precipitation associated with the HTM may enhance CAR, while wetter conditions associated with Neoglacial climate (post-3 cal. kyr BP) may suppress CAR. As a consequence of generally cooler, wetter conditions in the HBL during the mid- to late Holocene, we further suggest that such conditions may have supported long-term CH4 production through sustained minerotrophic to oligotrophic peatland conditions, a growing peat mass, and the potential for greater release of the products of anaerobic decomposition (Packalen et al., 2014). Accordingly, we examine peatland development and C dynamics in the HBL in relation to mid- to late Holocene atmospheric CH4 trends, and evaluate the potential for HBL peatlands to function as a late-Holocene terrestrial source of CH4.
Materials and methods
Study setting
The HBL is the second largest continuous peatland region globally (Riley, 2011), and is located between 50 N and 60 N, and 76 W, and 100 W (Figure 1). The northern and eastern boundaries of the HBL stretch to the margins of Hudson and James Bays, while the southern and western margins of the HBL follow the contour of the Canadian Shield. For the purposes of this study, the area occupied by the HBL is defined according to the Hudson Plains Terrestrial Ecoregion in the North American Environmental Atlas (http://www.cec.org/naatlas/). Using this definition, the HBL has a surface area of 372,000 km2, of which up to 90% is classified as peatland (Ontario Ministry of Natural Resources, 2000). The HBL landscape was shaped by the Laurentide Ice Sheet (LIS). Evidence of past glacial activity is present in the form of continuous and low-relief deposits of till and diamicton, eskers and moraines, and glacio-marine sediments that were deposited over top of Paleozoic sedimentary bedrock (Martini, 2006). Deglaciation resulted in high rates of GIA and subsequent flooding by the Tyrrell Sea which delayed land emergence and the onset of peat initiation until ~8 cal. kyr BP (Packalen et al., 2014). Furthermore, the final collapse of the LIS in the vicinity of the HBL also delayed the onset of the HTM for northeastern Canada until the mid-Holocene (Renssen et al., 2009).

Physical features and study locations within the physiographic region of the Hudson Bay Lowlands, Canada. The region includes the coastal Hudson Bay Lowland (HBL), inland HBL, and James Bay Lowland (JBL) ecoregions, where the boundary between the two regions is indicated by a thin black line (http://www.cec.org/naatlas/). Numbered contours (thick black lines) represent modern isostatic uplift rates (mm/yr; Peltier, 2004), while gray shading represents modern surface elevation generated using 1:250K Canadian digital elevation data (http://www.geobase.ca). Available basal dated peat core locations for estimating peatland expansion (n = 100; black dots; Packalen et al., 2014; and references therein) and vertically dated peat study locations for carbon analysis (n = 17, white triangles) are indicated.
Climate in the HBL is microthermal, strongly influenced by Arctic air masses and strong winds. Modern average climate (1971–2000) is characterized by an annual temperature of −2.5 ± 1.8°C (mean ± standard deviation (SD)), and total annual precipitation ranging from 429 to 743 mm. The growing season in the HBL stretches over a period of 119–162 days, and is characterized by a growing season temperature of 10.8 ± 1.0°C (mean ± SD), and summer precipitation, as rain, equivalent to half to two-thirds of the total annual precipitation (McKenney et al., 2006). Occasional sporadic to discontinuous permafrost features occur in the HBL peatlands, together with a narrow stretch of continuous permafrost limited to the northernmost reaches of the region, near the Hudson Bay coast (Riley, 2011).
Sample collection and data sources
During the 2009–2011 field seasons, complete peat profiles were collected from representative peatlands distributed across the HBL. To minimize peat compaction during sampling, surface peat (0–80 cm) was collected using a Wardenaar or Jeglum-type box corer (10 × 10 × 80 cm3), while the subsequent sampling to mineral contact was carried out using a Russian pattern side-cutting peat sampler (5 × 50 cm2). Recovered core segments were wrapped in plastic and aluminum foil and stored in polyvinylchloride pipe at −10°C until laboratory analysis.
For the purpose of examining C flux histories using the ‘Super Peatland’ approach described by Yu (2011), CAR was quantified using measurements of bulk density, C content, and vertical radiocarbon dating collected from 17 complete peat profiles. In all, 10 of these records are new contributions (14C-accelerator mass spectrometry (AMS)), while the remaining seven records were obtained from previously published literature and include both 14C-AMS and conventional dates (Bunbury et al., 2012; Kettles et al., 2000; Kuhry, 1998, 2008; O’Reilly et al., 2014). Peatland area increase in the HBL was estimated using a previously published synthesis of peat initiation history of the HBL, which comprises 100 new and previously reported peat basal ages (Packalen et al., 2014; and references therein). Data obtained from previously published literature were selected on the basis of three criteria: (1) their specified location within peatlands of the HBL ecoregion, (2) the availability of radiocarbon dates at the basal peat–mineral interface, and (3) the availability of total peat depth.
Peat physical properties
Bulk density, loss-on-ignition (LOI), and/or C content were obtained using standard methods (Yu, 2012, and references therein). Briefly, contiguous samples of known volume (5–10 cm3) were cut from fresh material at 2- to 4-cm intervals and oven-dried to constant mass. Direct measurement of C content was completed using an Elementar Variomax CN analyzer (Elementar Analysensysteme GmbH Donaustraße 7 63452 Hanau Germany). Data from the seven previously reported peat records were re-interpreted using author-contributed raw datasets of bulk density, LOI and directly measured C content (Bunbury et al., 2012; O’Reilly et al., 2014), or C content at a rate of 50% LOI, where measured C content was not available (Kettles et al., 2000; Kuhry, 1998, 2008).
Bulk density (g/cm3) was calculated by dividing the dry peat mass (g) by the fresh peat volume (cm3), while C density (g/cm3) was calculated by multiplying the bulk density of each peat increment by the corresponding C concentration. Areal dry peat mass (g/cm2) was determined as a function of peat depth, and cumulative dry mass was determined on an area basis as, averages weighted by layer thickness. Here, we consider the rate of C accumulation using several conventions commonly reported in the literature, including (1) the long-term apparent rate of C accumulation (LORCA), representing average C accumulation since peat inception; (2) CAR, which uses vertical chronological control to track temporal patterns in the rate of C accumulation; and (3) NCB, the peat decay adjusted product of peatland area and CAR (discussed in detail below), which models the net C sequestration by peatlands through time and is potentially the most appropriate term to use when comparing historical peatland C dynamics and the global C cycle (Yu, 2012).
Radiocarbon dating
New radiocarbon dates were obtained using identified above-ground plant macrofossils sampled from a 1-cm peat segment from 10 cores (Table 1). 14C-AMS dating of these samples was completed at either the Keck-CCAMS facility (Irvine, USA), the UGAMS facility (Athens, USA), or at Beta Analytic, Inc (Miami, USA). All dates compiled from previous studies, which include 14C-AMS dated peat (Bunbury et al., 2012; O’Reilly et al., 2014), and both 14C-AMS and radiometric dates (Kettles et al., 2000; Kuhry, 1998, 2008), in addition to the newly dated samples, were calibrated using IntCal09 (Reimer et al., 2009) in the Clam package for R (Blaauw, 2010). All calibrated ages are expressed as calendar years before present (cal. yr BP), where present is
14C-accelerator mass spectrometry (AMS) dating of peat macrofossil for 10 new sites in the Hudson Bay Lowlands, Canada; by increasing latitude.
Calibrated with IntCal09 calibration curve (Reimer et al., 2009).
Peat decay
To estimate the true rate of peat C accumulation, the long-term rate of C decay in the peat profile must be quantified. Other studies where this model has been applied have concluded that the constant decay function provides the simplest best fit to the data (Yu, 2011). Accordingly, the exponential bog decay model (Eq. 1) was applied to the subset of 17 HBL peat records (Clymo et al., 1998). This model states that:
where M is the cumulative peat C mass (g C/cm2), p is the annual peat addition to the catotelm and determines the slope of M versus peat age curve, a is the peat decay constant that determines the curvature of M versus peat age, and t is time. The subset of new and/or synthesized age–depth models was fitted using Eq. 1, such that depth was represented as cumulative peat C mass (g/cm2). Alternative and potentially more ecologically meaningful peat decay rules, such as linear and non-linear (quadratic) decay through time were also explored (Clymo et al., 1998). Although they appeared to impact estimates of total C losses, they were not found to not improve curve fitting and were not used in subsequent analyses.
Modeling of peatland C dynamics
Peatland C dynamics since mid-Holocene peat initiation for the HBL was quantified following the methods of Yu (2011), using mean CAR from the 17 peat profiles described above, binned at 0.5-kyr intervals. Table 2 provides a summary of the terms used to describe peatland C dynamics in the HBL; however, additional details are presented by Yu (2011). Briefly, NCB refers to the difference between NCU at a given time interval and net C release (NCR) by the whole peatland, due to decay, prior to the given time interval. To quantify NCB, mean time-weighted CAR (n = 17) were determined for 0.5-kyr peat cohorts. The net C pool (NCP) was then quantified for each peat cohort by multiplying the mean CAR at each 0.5-kyr bin by the corresponding total peatland area at that time. Peatland area at each time interval was estimated as a relative proportion of cumulative peat initiation, assuming a modern HBL extent of 372,000 km2, with 90% occupied by peatlands. NCU for each peat cohort was calculated as the total amount of C at each time interval, accounting for previous peat losses because of decomposition using the exponential decay function calculated above (NCP*1/e−at). Finally, NCR was quantified as the total amount of C released by the entire peatland, including previously accumulated peat, prior to the relative time period. Cumulative C pool, uptake, release, and balance were also calculated as cumulative sums of the respective net C terms to track the total C terms through time in the HBL.
Summary of terms used to describe peatland carbon (C) dynamics in the Hudson Bay Lowlands, Canada (adapted from Yu, 2011).
All terms reported in Pg C.
Results
C accumulation and peat decay rates
Vertical radiocarbon data for 10 new peat profiles are presented in Table 1, and were combined with an additional seven vertically dated peat profiles drawn from the published literature (Bunbury et al., 2012; Kettles et al., 2000; Kuhry, 1998, 2008; O’Reilly et al., 2014). With the exception of two coastal peatland records from the north-west HBL, a region of permafrost influence, most records considered here are located in the inland peatland region and include bog and fen peatlands. Shallow coastal peatlands were not sufficiently well constrained in terms of their chronologies for inclusion in the present analysis. With that in mind, the total peat depth of the 17 vertically dated peat records ranged from 131 to 311 cm, while basal peat ages ranged between 3940 and 6810 cal. yr BP.
C content per unit dry peat varied along the peat profile from 12% to 54%, with the lowest values being recorded near the mineral contact. The mean and standard error (SE) of peat C content (n = 885 from 17 cores) was 48 ± 0.2%, and bulk density (n = 885 from 17 cores) was 97.4 ± 2.3 g/dm3 (n = 885 from 17 cores). Using only the basal age of the peat deposit for the subset of peat cores used in this study (n = 17), we found that median LORCA ranged between 11.2 and 28.3 g C/m2/yr (median: 18.3 g C/m2/yr). Consideration of temporal variability in peat accumulation rates yielded CAR estimates that varied over two orders of magnitude among intervals, and ranged from 3.5 to 310 g C/m2/yr. Incorporating vertical radiocarbon peat dates results in a Holocene CAR of 25 ± 0.9 g C/m2/yr (SE; median: 18.5 g C/m2/yr). Upon closer inspection of the peat C accumulation records, temporal variability in CAR appears to be related to periodic increases in age–depth model slopes through time, coupled with greater C content in fen peatland records, which has been observed elsewhere (Holmquist and MacDonald, 2014; Loisel et al., 2014).
Peat accumulation in the HBL varied among the 17 sites examined in this study, with concave, convex, and linear patterns of peat age versus peat accumulation observed (n = 80 from 17 cores). Considering all sites together, we found a linear relationship (Figure 2a) between peat depth and cumulative peat OM mass (R2 = 0.74, p < 0.0001, n = 80), as well as between peat depth and cumulative peat C mass (R2 = 0.93, p < 0.0001, n = 80). Fitting the C mass data with the simplest decay model (Eq. 1), which assumes constant decay through time, to the relationship between peat age and cumulative C mass for all HBL sites together (n = 80 from 17 cores) yielded a modeled peat C addition rate (PCAR; Yu, 2011) of 21 g C/m2/yr and a decay constant, a = 0.0000562/yr for the HBL (Figure 2b). Though not statistically significant in this study, application of alternative peat decay rules (Clymo et al., 1998), such as linear and non-linear (quadratic) decay through time increased the decay coefficients to 0.000067/yr and 0.000082/yr, respectively. Thus, for the purposes of this study, the simplest decay term was selected; however, further investigation is warranted into ecologically meaningful representations of decay in peatlands of differing trophic states, the results of which may better constrain total C mass estimates in the HBL.

Peat age, depth, and cumulative carbon (C) mass relationships and modeled exponential peat decay for the patterned peatlands of the Hudson Bay Lowlands (HBL), Canada. (a) Linear relationships among median peat age (n = 80 from 17 cores; dots) and 2σ age ranges (error bars) are presented for 17 vertically dated peat cores sampled from the HBL. Peat ages were estimated following calibration of 14C ages using the IntCal09 curve (Reimer et al., 2009), and (b) cumulative C mass was fit with an exponential peat decay model (Clymo et al., 1998). The best-fit peat C addition rate (PCAR) and long-term peat decay parameter (a) for the HBL are presented.
Modeled peatland C dynamics in the HBL
NCB (Figure 3a) is a spatially constrained snapshot of net C sequestration by representative peatlands in the HBL using the most comprehensive regional record of peat initiation and expansion currently available. The sum of modeled NCB (CCB) must be equal to the sum of the measured NCP (CCP; Yu, 2011); however, the trajectory in arriving at the total cumulative C term for each differs since peat initiation, owing to the accounting of delayed peat decay in the NCB term. According to the assumptions and limitations described above, the modeled potential NCU by HBL peatlands during the mid- to late Holocene ranged from 0.1 to 6.5 Pg C per 0.5-kyr interval, with a cumulative C uptake (CCU) since peat initiation of 45 Pg C; assuming constant peat decay since peat inception and across the growing peat landscape revealed a potential loss (CCR) of 7.6 Pg C (Figure 3b). The difference between NCU and NCR produces an NCB ranging from 0.1–5 Pg C, which summed results in a modeled total HBL C pool of 37 Pg C. Consideration of alternative decay rules, such as non-linear peat decay, which may be more ecologically meaningful, resulted in the largest decay coefficient for the HBL (a = 0.000082/yr, presented above), and nearly doubles the total potential peat C loss to ~15 Pg C. Application of the non-linear decay coefficient to NCU reduces the modeled NCB for the HBL to ~30 Pg C, more consistent with previous estimates using LORCA and peatland expansion dynamics to estimate the magnitude of the HBL C pool (Packalen et al., 2014). The outcome of curve fitting using non-constant decay terms yielded non-significant co-variates, thus such results are presented here for qualitative purposes alone to highlight the need for an improved understanding of peat decay dynamics in the HBL.

Holocene peat carbon (C) pools and modeled peat C terms (Yu, 2011) for the Hudson Bay Lowlands, Canada, since peat initiation began ~8 cal. kyr BP. (a) Peat age versus measured (n = 17 cores) net peat C pool (NCP) and modeled (n = 100 cores) net C uptake (NCU), release (NCR), and balance (NCB), where NCB = NCU − NCR., and (b) peat age versus cumulative measured peat C pool (CCP) and modeled cumulative peat C uptake (CCU) and release (CCR) following long-term constant peat decay, are presented at 0.5-kyr intervals.
The peatlands of the HBL have been a persistent C sink since they began to develop during the mid-Holocene (Figure 4). During the past 7.5 kyr of peat initiation and expansion in the HBL, mean CAR (0.5-kyr bins) ranged from 14 to 38 g C/m2/yr (Figure 4d). Late Holocene CAR in the HBL remained relatively stable, below the mean CAR (25 g C/m2/yr). By contrast, mid-Holocene CAR was highly variable in the HBL, and peaked (~35 g C/m2/yr) twice during this period: first 6–6.5 cal. kyr BP, and again 4–4.5 cal. kyr BP, following a 1.5 cal. kyr period of below average peat CAR. The variation in CAR appears to track patterns in peatland initiation, such that maximum CAR occurred synchronously with peat initiation maxima (Figure 4e). CAR also peaks in the most recent past (0.5 cal. kyr BP), likely related to under-decomposed acrotelm peat. Regional pollen-inferred climate records reveal relatively cold and dry conditions versus contemporary climate and a trend toward warmer and wetter conditions between 7.5 and 4.5 cal. kyr BP, when peat CAR was greatest and most variable. Mid-Holocene CAR maxima also appear to be preceded by warmer summers and wetter conditions (Figure 4c), which may suggest that periodic enhanced peatland productivity may contribute to CAR variability. By contrast, relatively similar to slightly cooler conditions compared with modern climate coupled with anomalously wet conditions were associated with the lowest and least variable CAR, though some minor variation is noted.

Holocene peatland area increase, carbon (C) accumulation, and net C balance, since peat inception for the Hudson Bay Lowlands (HBL), Canada. Data are presented in relation to regional paleoclimate, global insolation trends, and atmospheric methane (CH4) and carbon dioxide (CO2) variability. (a) Winter and summer insolation during last 7.5 kyr at 60°N (Berger and Loutre, 1991); (b) GISP2 (Greenland) ice core record of atmospheric CH4 concentration (Brook et al., 2000) and Antarctic ice core record of atmospheric CO2 concentrations (Monnin et al., 2004); (c) reconstruction of winter and summer temperature and annual precipitation anomalies inferred using a synthesis of the North American Pollen Database for the northern Quebec, Canada pollen region (Viau and Gajewski, 2009); (d). Mean (± standard error) time-weighted C accumulation rates (CARs, n = 17), and modeled (Yu, 2011) net carbon uptake and balance, at 0.5-kyr intervals (n = 100); (e) frequency of post-glacial peatland initiation in the HBL, inferred from basal radiocarbon dates (n = 100, 0.5-kyr bins), cumulative peatland initiation (%; n = 100), and peatland surface area increase (%, n = 100; Packalen et al., 2014).
We calculated a modeled peat C sequestration rate for the HBL of ~7 Tg C/yr and a mean NCB of ~4.4 Tg/yr. Instantaneous CARs generated from the NCB term account for both C uptake and C release over the period of peatland development (Yu, 2012). In the HBL, the modeled instantaneous CARs ranged between 14 and 42 g C/m2/yr, and were greatest during the mid-Holocene, coincident with the period of rapid land emergence from the Tyrrell Sea and associated intense peat initiation. A delay in peat decay is apparent in the HBL, such that most of the modeled C was released during the late Holocene. Accordingly, our results indicate that ~6.4 Pg C was released between 0 and 4 cal. kyr BP, compared with a modeled mid-Holocene (4–7.5 cal. kyr BP) release of ~1.2 Pg C, for a total of 7.6 Pg C since the time of peat initiation in the HBL.
Discussion
Carbon accumulation patterns in the HBL
Peat accumulation in the HBL has occurred on a topographically favorable landscape since the mid-Holocene, against a backdrop of decreasing summer insolation and insolation seasonality, and within climatic boundary conditions suitable for peatland development (Beilman et al., 2009; Charman et al., 2013; Packalen et al., 2014; Viau and Gajewski, 2009). GIA-driven land emergence in the HBL, which decays exponentially through time, is the fundamental control on the timing of peat initiation in this region (Packalen et al., 2014). LORCA is a widely reported metric for comparing mean C accumulation among peatland regions (Tolonen and Turunen, 1996). Previous results estimated using a synthesis of 100 basal peat ages reveal a mean ± SD LORCA equal to 18.5 ± 5.7 g C/m2/yr for the HBL (Packalen et al., 2014). This rate is within the range of LORCA found for the subset of peat cores used in the present study (median: 18.3 g C/m2/yr; n = 17).
LORCA in the HBL is comparable to estimates reported for similarly aged peatlands, such as the Eastmain region, Quebec, Canada, where inundation by the Tyrrell Sea delayed peat initiation until the mid-Holocene, and C accumulated from 16.2 to 18.5 g C/m2/yr (Van Bellen et al., 2011). Furthermore, the HBL LORCA is also similar to the rates reported for older peatland regions that began to initiate in the early Holocene, such as the West Siberian peatlands where the rates range from 12.1 to 23.7 g C/m2/yr (Turunen et al., 2001), undrained Finnish peatlands with a rate of 18.5 g/m2/yr (Turunen et al., 2002), and the continental peatlands of western Canada with a rate of 19.4 g C/m2/yr (Vitt et al., 2000).
The rate of peat C accumulation following peat inception is related to both autogenic and allogenic factors, and likely varies through time. However, LORCA as it was originally defined (Tolonen and Turunen, 1996) may not adequately account for peat mass losses through time due to fire and/or erosion, temporal gaps in C accumulation due to permafrost accretion or enhanced decomposition, and other disturbances (Clymo et al., 1998; Tarnocai et al., 2012; Yu, 2012). Consequently, LORCA does not permit a complete assessment of temporal patterns in C accumulation in relation to changes in peat burial rates related to peat type and oligotrophication, growing season length, and/or other potential environmental controls (Loisel et al., 2014).
As an alternative, we also consider CAR, a metric that captures the apparent temporal variability in the rate of C accumulation, through vertical chronological control using multiple radiocarbon dates along the peat column. Using the subset of well-dated cores from the HBL (n = 17), our results reveal a Holocene CAR of ~25 g C/m2/yr (median: 18.5 g C/m2/yr) for the HBL, which is similar to the recently reported northern CAR of 22.9 ± 2 g C/m2/yr (SE; Loisel et al., 2014), and to values from nearby subarctic Quebec, Canada (~24 g C/m2/yr; Lamarre et al., 2012), and neighboring boreal shield peatland region of northwestern Ontario, Canada (~24 g C/m2/yr; Holmquist and MacDonald, 2014).
Temporal variability in CAR may be related in part to the vegetation community composition and hence relative decomposability of the plant material through time. Previously, bog peat dominated by Sphagnum species was shown to be less C-rich than peat associated with more minerotrophic conditions (Beilman et al., 2009). Furthermore, periods of increased CAR identified for the northern peatland region (Loisel et al., 2014), were attributed to changes in plant composition and C density through time, rather than increasing peat bulk density because of compaction of older peat cohorts. In the HBL, plant macrofossil and pollen-based peat vegetation reconstructions provide evidence that minerotrophic conditions following peat inception can persist for millennia (Bunbury et al., 2012; Holmquist and MacDonald, 2014; O’Reilly et al., 2014; Sjörs, 1959; Van Bellen et al., 2013). Thus, for the 17 peatlands examined here, C-rich inputs from young, minerotrophic peatland vegetation, coupled with periods of rapid rates of peat accumulation, may contribute to explaining apparent CAR fluctuations in this record.
Peat initiation frequency in the HBL was most intense during the mid-Holocene followed by a reduction in the rate of peat initiation through the late Holocene (Packalen et al., 2014). This peat initiation pattern of generally declining rates though time is associated with the GIA-controlled reduction in new land available for occupation by new peatlands. Comparing these trends with HBL CAR reveals that the mid-Holocene periods of intensified peat initiation in the HBL were accompanied by elevated CAR of ~35 g C/m2/yr, with maxima around 4–6 cal. kyr BP (Figure 4). As peat initiation rates declined through the late Holocene, partly in response to exponentially declining rates of GIA, late Holocene CAR was weakly variable and remained below the long-term mean CAR.
By the onset of the late Holocene (4 cal. kyr BP), ~75% of the HBL peatlands had initiated, and some of these peatlands were on a trajectory toward lower C density, Sphagnum-dominated, oligotrophic conditions. Previous evidence, inferred using testate amoebae and plant macrofossils, indicates that many modern bogs in the vicinity of the Hudson Bay transitioned from the fen stage in the last millennium or during latter half of the late Holocene (Holmquist and MacDonald, 2014; Lamarre et al., 2012). Today, ~40% of the HBL remains classified as fen peatlands (Ontario Ministry of Natural Resources, 2000). Thus, lower rates of peat initiation in the HBL during the late Holocene, and the corresponding trend of lower CAR during this period may be more strongly influenced by peat succession dynamics occurring in existing HBL peatlands, than by C inputs from newly emerging peatlands.
Modeled peatland C dynamics in the HBL
Using all available peatland records for the HBL with vertical chronological control, which included minerotrophic and oligotrophic peatlands in both the presence and absence of permafrost, we assess the temporal variation in Holocene C uptake and release using the ‘Super Peatland’ approach (Yu, 2011). As part of this process, peat C decay was estimated assuming constant decay and resulted in a slightly concave decay model of age versus depth. Though frequently considered, constant decay may be an oversimplified assumption, as individual peatlands or peat types may exhibit very different decay patterns or may not be well fitted to the decay model described by Clymo et al. (1998). In addition to this limitation, young, shallow coastal peatlands are not well represented in this study, owing to the challenges of dating the relatively young peat deposits. Thus, we also consider alternative decay assumptions qualitatively, as model outputs yielded non-significant terms in the present study, and sample limitations precluded additional investigation. Nevertheless, using a constant peat decay term (Clymo et al., 1998) summation of the NCB model output term, resulted in a total contemporary C mass equal to 37 Pg C (Figure 4). Under conditions of constant peat decay, this modeled total C mass estimate is ~20% greater than a previous estimate of 30 Pg C obtained using patterns of peat initiation and expansion, land emergence, and measured C mass (Packalen et al., 2014). Still, peat decay is not well understood for HBL peatlands, especially regarding the relative role of fen versus bog peatlands. Consideration of alternative decay rules in the calculation of NCB, such as non-linear peat decay, nearly doubles the total C release term (NCR), and results in a reduced HBL C mass estimate of 30 Pg C, consistent with previous estimates (Packalen et al., 2014). However, we note here that the non-constant decay terms were not significant in this study, and thus are used only for qualitative purposes to demonstrate the effect of enhanced peat decay, at the landscape scale, on the total C pool in the HBL. Our results confirm that a better understanding of peat decay, especially concerning the influence of fen versus bog peat cover, is needed to fully explore Holocene C dynamics.
Although the absolute values for HBL C-pools reported here are uncertain, temporal patterns in Holocene C dynamics may be more meaningful and suggest that a total of 7.6 Pg C may been lost from the total peatland C mass since peat inception began in the HBL during the mid-Holocene. Maximum modeled C sequestration and release was greatest during the late Holocene, as the peat complex expanded across the rapidly emerging low-relief landscape. Additional decomposition of C flux terms reveals a mid-Holocene (4–7.5 cal. kyr BP) net HBL C sequestration equal to 14 Pg C, which occurred at a mean rate of 4 Tg C/yr, and included a small C release of 1.2 Pg C. During the late Holocene, new peatland initiation coupled with rapid peat accretion in existing peatlands resulted in an additional C sequestration of 23 Pg C, at a mean rate of 6 Tg C/yr, and includes a larger potential C release of 6.4 Pg C. The larger late Holocene C release accounts for the gradual decay in the catotelm of previously deposited peat cohorts, which is expected (Yu, 2011), although the effect of decay is smaller in the HBL, due in part to the relatively young age of the HBL peat complex.
NCB and paleoclimate in the HBL
The timing of peat initiation events in the HBL is primarily driven by uplift; however, paleoclimate may have contributed to facilitating initiation events, and perhaps even more so, in influencing the rates of Holocene C accumulation and NCB in the millennia following initiation. As a consequence of the delayed retreat of the LIS near the HBL, the timing of the HTM in northeastern Canada is thought to have been deferred until ~6–7 cal. kyr BP (Renssen et al., 2009). The delay of the HTM in the HBL until the mid-Holocene corresponds to the earliest period of intensified peat initiation and a period of peak C accumulation in the HBL. Coupled climate model simulations and pollen reconstructions for north-east Canada suggest a possible warming of ~1°C versus pre-industrial temperatures, which may account for higher C sequestration rates via enhanced primary production during this time. However, regional paleoclimate reconstructions appear to confound this relationship.
Comparisons using apparent CAR among Canadian peatlands regions in the vicinity of the HBL reveal that maximum CAR occurred during the mid-Holocene, while the lowest CAR occurred in the late Holocene (Van Bellen et al., 2011; Vitt et al., 2000), and early Holocene for the continental peatlands of western Canada. Lower CAR in both the western Canadian peatlands and the Eastmain region, Quebec, during the late Holocene was attributed to a cooling climate and drier peatland conditions related to height-induced surface drying and/or permafrost accretion (Van Bellen et al., 2011; Vitt et al., 2000), while low early Holocene CAR was attributed to dry peatland conditions (Vitt et al., 2000). However, pollen-inferred climate reconstructions for northern Quebec, Canada (Viau and Gajewski, 2009), that are consistent with an independent paleoclimatic record from the HBL (McAndrews et al., 1982), imply lower than present temperature and precipitation during the mid-Holocene (until 4.5 cal. kyr BP). Thus, in the HBL, the negative precipitation anomaly (i.e. reduced precipitation compared with the contemporary climate normal) corresponded to higher CAR during the mid-Holocene. Given the uncertainty in mid-Holocene climate in the HBL, the relationship between CAR and climate is interpreted cautiously here and warrants further investigation.
Regional paleoclimate reconstructions for the late Holocene near the HBL reveal climatic conditions similar to today, with the exception of the last 1.5 kyr when increased precipitation and reduced temperatures associated with Neoglacial cooling are recorded. During this period, CAR appears to increase slightly as precipitation increases, reaching a late Holocene peak CAR around 1.2 cal. kyr BP followed by a small decrease in CAR around 0.8 cal. kyr BP, around the ‘Medieval Climate Anomaly’ (MCA). Similarly, a decline in CAR during the MCA is also apparent in the northwestern Ontario peatlands, adjacent to the HBL (Holmquist and MacDonald, 2014), while a slight increase is noted at the time of the ‘Little Ice Age’ (LIA). Evidence recorded in Quebec, Canada, peatlands to the east of Hudson Bay document a wet shift during the LIA and a return of minerotrophic conditions, which were attributed to a cooler, wetter climate (Van Bellen et al., 2013). However, our conclusions in this regard are highly limited by the resolution of the time series used for CAR in the present study (0.5-kyr bins) compared with the century scale of the MCA and LIA climatic events. This limitation highlights the need for further high resolution examinations of these time periods to better constrain the relationship between climate and CAR in the HBL.
Contemporary surface measurements have revealed that C losses from peatlands are tightly coupled with vegetation composition, trophic status, and hydroclimatic conditions (Bubier, 1995; Moore and Dalva, 1993), and can occur via aerobic decomposition in the acrotelm and by anaerobic decomposition in the catotelm (Clymo et al., 1998; Gorham, 1991). Furthermore, rapid initiation and expansion of minerotrophic fens in the early Holocene have been proposed to have contributed to peak early Holocene CH4 concentrations (MacDonald et al., 2006; Smith et al., 2004). Further, sources of early Holocene CH4 emissions have more recently been expanded to include pre-peatland wetlands such as marshes and wet fens as facilitators of biosphere CH4 transport to the atmosphere (Yu et al., 2013). In the HBL, surface and airborne measurements of contemporary snow-free CH4 flux, interpreted using a chemical transport model (Pickett-Heaps et al., 2011) suggest that the largest rates of net CH4 flux occur in association with the older, inland peatland region rather than the younger coastal region. Accordingly, Packalen et al. (2014) recently hypothesized that the largest atmospheric contributions of peat-derived CH4 from the HBL may have occurred during the late Holocene, and suggested that a 3 cal. kyr BP peat mass could potentially release 1–7 Tg CH4/yr to the late Holocene atmosphere.
Here, we examined patterns in Holocene NCB of HBL peatlands, for comparison with ice core inferred atmospheric trace gas variation (Figure 4). Similar climate-CAR patterns are also apparent in northwestern Ontario peatlands, adjacent to the HBL (Holmquist and MacDonald, 2014). Assuming a constant decay term, our model estimations predict a late Holocene C release of 6.4 Pg C, equivalent to a mean rate of 1.6 Tg C/yr (range: 1–2.5 Tg C/yr) from the more developed HBL peatland landscape. Alternative non-linear decay patterns, though not significant in this study, but may be more ecologically realistic, suggest even larger C losses from the HBL C mass, such that total late Holocene C losses increase to ~10.2 Pg C, corresponding to a mean rate of ~2.6 Tg C/yr (range: 1.5–4.1 Tg C/yr) during the late Holocene and from across the HBL. Assuming a portion of these C losses occur as CH4, then this model evidence, irrespective of decay assumption is comparable to previous estimates of late Holocene CH4 production potential from the HBL. Therefore, our analysis of modeled C dynamics provides an additional line of evidence of natural terrestrial C contributions from the HBL to the late Holocene atmosphere, particularly in the form of CH4 flux, as a consequence of the minerotrophic to oligotrophic peat patterning that characterizes the HBL landscape.
The trajectory of NCB in the HBL differs from other important peatland regions globally, such that C sequestration was greatest during the late Holocene, rather than at a minimum. While evidence of peat succession is present in the HBL, as peatlands shift from minerotrophic to oligotrophic systems, our data show that the relationship between cumulative peat mass and peat depth in the HBL remains linear. Assuming that the relationship between peat mass and depth is related in part to trophic succession in the HBL, and given that there is little evidence to suggest that mass peat losses are occurring as a consequence of drying and oxidation of peat, the peatlands of the HBL may continue to accumulate peat (and sequester C) for some time to come. Our findings convey the importance of the HBL in global C accounting, and suggest that significant C losses from the HBL may have important implications for the global C budget and climate system.
Conclusion
Here, we have examined trends in the timing and magnitude of variation among CAR and net C uptake, release, and balance since peat inception began in the HBL. Previous findings have confirmed the fundamental control of GIA over the timing of peat initiation in the HBL and that peatland development has occurred within the context of a weakly varying climate (Glaser et al., 2004a, 2004b; Packalen et al., 2014). Our findings do not reveal a tight coupling between climate and CAR in the HBL; however, they do suggest a potential link with climate that warrants further investigation. Furthermore, we report that CAR may be more strongly related to peatland succession dynamics in the HBL. Accordingly, we find that CAR is related to the intensity of peat initiation, the latter of which determines the proportion of young peatlands with high values for CAR on the rapidly emerging HBL landscape. Our data show that, while the HBL has been a persistent C sink for millennia, more than two-thirds of the total C mass accrued during the late Holocene, when CAR remained below the long-term mean. Most of the HBL peatlands initiated during the mid-Holocene, in advance of the late Holocene atmospheric CH4 rise. Moreover, our findings suggest that most of the potential C lost from the HBL occurred during the late Holocene, likely owing to decay of previously deposited peat. Modern evidence suggests the largest CH4 emissions in the HBL occur in association with older, inland patterned peatlands (Pickett-Heaps et al., 2011). Given that sustained minerotrophic to oligotrophic peat patterning typifies the HBL landscape, the persistent C release during the late Holocene may provide a first line of evidence of natural terrestrial C contributions from the HBL to the late Holocene atmosphere, particularly in the form of CH4.
Footnotes
Acknowledgements
Sincere thanks to Jim McLaughlin for research funding and field support. We also thank A Dyke for providing access to the Canadian basal radiocarbon database and P Kuhry for contributing raw peat core data.
Funding
The study was partially funded by Ontario Ministry of Natural Resources’ Applied Research and Development Branch and Far North Branch, under the auspices of projects CC-021 and FNIKM 028. Additional support for field work and radiocarbon dating was provided by grants from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Ontario Ministry of the Environment through the Climate Change and Multiple Stressor Research Program at Laurentian University to SAF. Support in the form of an NSERC Alexander Graham Bell Canada Postgraduate Scholarship and an Association of Canadian Universities for Northern Studies, Canadian Northern Studies Trust Scholarship, and field research grants from the Society of Wetland Scientists and Aboriginal Affairs and Northern Development Canada’s Northern Scientific Training Program were provided to MSP.
