Climate-related changes in peatland carbon accumulation during the last millennium

. Peatlands are a major terrestrial carbon store and a persistent natural carbon sink during the Holocene, but there is considerable uncertainty over the fate of peatland carbon in a changing climate. It is generally assumed that higher temperatures will increase peat decay, causing a positive feedback to climate warming and contributing to the global positive carbon cycle feedback. Here we use a new extensive database of peat proﬁles across northern high latitudes to ex-amine spatial and temporal patterns of carbon accumulation over the past millennium. Opposite to expectations, our results indicate a small negative carbon cycle feedback from past changes in the long-term accumulation rates of northern peatlands. Total carbon accumulated over the last 1000 yr is linearly related to contemporary growing season length and photosynthetically active radiation, suggesting that variability in net primary productivity is more important than decomposition in determining long-term carbon accumulation. Furthermore, northern peatland carbon sequestration rate de-clined over the climate transition from the Medieval Climate Anomaly (MCA) to the Little Ice Age (LIA), probably because of lower LIA temperatures combined with increased cloudiness suppressing net primary productivity. Other factors including changing moisture status, peatland distribution, ﬁre, nitrogen deposition, permafrost thaw and methane emissions will also inﬂuence future peatland carbon cycle feedbacks, but our data suggest that the carbon sequestration rate could over many areas of northern peatlands in a warmer


Introduction
Peatlands contain around 600 gigatonnes of carbon (Gt C) that has accumulated since the last glacial maximum in northern mid-high latitudes, tropical regions and temperate areas of the Southern Hemisphere, and the steady accumulation of carbon has been a small but persistent sink for atmospheric CO 2 throughout the Holocene (Yu, 2011). The relationship between climate change and the rate of carbon sequestration is important for understanding the past and future global carbon cycle, and it has generally been assumed that because temperature drives increasing decay (Ise et al., 2008;Dorrepaal et al., 2009), peatlands could be part of the positive feedback from the global carbon cycle (Friedlingstein et al., 2006). A key objective in improving understanding of the global carbon cycle in climate models is to be able to simulate past observed atmospheric CO 2 changes.
There is growing interest in the last millennium as a climate-modelling target, and especially in the assessment of the sensitivity of the global carbon cycle to climate warming (Abe-Ouchi and Harrison, 2009;Jungclaus et al., 2010). In the Northern Hemisphere, the transition from the generally warmer Medieval Climate Anomaly (MCA) to the cooler Little Ice Age (LIA) (Mann et al., 2008(Mann et al., , 2009Jansen et al., 2007) was associated with a ca. 7-10 ppmv decline in atmospheric CO 2 concentration (Ahn et al., 2012). This pattern supports the existence of a positive global climate-carbon cycle feedback, as suggested by coupled climate-carbon cycle models (Friedlingstein et al., 2006;Denman et al., 2007). However, estimates of the magnitude of the climate sensitivity of the global carbon cycle based on data from the last millennium vary from 1.7-21.4 ppm CO 2 K −1 (Frank et al., 2010) to 40-60 ppm CO 2 K −1 (Cox and Jones, 2008). Carbon cycle models also vary greatly in their assessment of this feedback (Friedlingstein et al., 2006), although recent estimates (Jungclaus et al., 2010) suggest sensitivity within the lower end of this range (3.2-12 ppm CO 2 K −1 ). The causes of the reduction in CO 2 concentrations during the MCA to LIA transition are poorly known, but reduced soil heterotrophic respiration is assumed to be important (Jungclaus et al., 2010;Pongratz et al., 2009). However, the models do not specifically take into account possible climate-related variations in the rate of peatland carbon sequestration.
Peatlands have sequestered and exchanged atmospheric carbon over millennia (MacDonald et al., 2006;Frolking and Roulet, 2007), with the largest store in northern extratropical peatlands, an estimated 545 Gt C (Yu et al., 2010). The annual uptake of CO 2 by peatlands, previously estimated as 0.076 Gt C yr −1 (Gorham, 1991) or 0.088 Gt C yr −1 (Yu, 2011), without considering long-term decay (see below), is a small but temporally persistent component of land carbon uptake. This is equivalent to 36 ppm atmospheric CO 2 over 1000 yr, based on a simple conversion from change in carbon pool to atmospheric CO 2 of 1 Gt C = 2.123 ppm. However, over millennial timescales, carbon uptake of this magnitude would be compensated by ocean outgassing processes as a result of reduced CO 2 in the atmosphere and reduction in air-sea CO 2 partial pressure, so that the actual effect on the atmosphere is only 20-35 % of this total over periods of 200-2000 yr (Archer et al., 2009), or 7-12 ppm atmospheric CO 2 over a 1000 yr period. Variations in the size of the peatland sink could therefore have a significant cumulative effect on global atmospheric CO 2 concentrations over the last millennium, of the same order of magnitude as the observed changes.
In this study, we compiled peat core data from northern peatlands to estimate changes in carbon accumulation over the last millennium and to explore the spatial relationship between climate and the total size of the carbon sink accumulated over this period. Further analysis on a subset of well-dated cores allowed an analysis of temporal variation in carbon accumulation in relations to the MCA-LIA climate changes estimated from palaeoclimate records. We use these data to help understand the relationships between climate and peatland carbon accumulation and to assess the direction and strength of the peatland carbon cycle feedback.

Site selection and carbon measurement
A list of Northern Hemisphere, extratropical peatland profiles with published and unpublished carbon accumulation data was compiled (Tables 1 and 2) for sites that met the following criteria: a. at least 3 evenly spaced dates (including 210 Pb, tephra, spheroidal carbonaceous particles, pollen markers or 14 C (pre-or postdating the period of nuclear bomb testing), and the uncut peat surface) and spanning the last ca. 1000 yr. Most of the sites had more than 5 dates (Table 2), but we also rejected some sites where a satisfactory age-depth model could not be produced, because of age reversals or other problems; and b. contiguous bulk density measurements at < 5 cm resolution.
Application of these criteria resulted in the selection of 24 sites that were used in subsequent analyses of the temporal changes in carbon accumulation through the last millennium (Table 2). A second tier of sites (Table 1) was used for the millennium carbon inventory analysis against climate indices. These sites did not meet criteria (a) and (b), but did meet the following criteria: c. a basic age-depth model for the last millennium; and d. contiguous bulk density measurements but not necessarily at high resolution.
A total of 90 sites met these less stringent criteria and were used in the climate-carbon inventory analyses (Tables 1 and 2). The sites in both data sets are widely distributed geographically and broadly representative of the climate space occupied by northern peatlands (Fig. 1).  Bulk density was measured on carefully cut fresh or frozen material using freeze drying or oven drying of samples of known volume. Sample sizes varied depending on the sampling method and core size, and sample resolution varied from 0.5 to 5 cm 3 (Table 2). In all cases samples were large enough to accurately measure bulk density and were taken contiguously to enable reliable estimates of dry mass accumulation over time. Carbon density was derived from bulk density multiplied by the carbon content for each sample. Where carbon data were not available, we assumed that 50 % of the organic fraction (measured by standard losson-ignition analysis at 500 • C) was organic carbon. A carbon value of approximately 50 % is routinely used for peat (Gorham, 1991;Vitt et al., 2000) and is reasonable compared to the mean carbon content of the nine sites for which we have measured values in this study (46.6 ± 0.33 %), and other studies in western Canada (51.8 %; Yu et al., 2009) and West Siberia (50.7-56.3 %;. To provide an assessment of hydrological differences among the peatlands in our analyses, we classified sites as either bogs or fens. Although differences between these two peatland types are related to the relative influence of different water sources (i.e. groundwater, surface water, precipitation), thresholds used for distinction between the two types are regionally varied. For our site classification, we used a relatively conservative approach, including only Sphagnumdominated systems that lacked vegetative or morphological evidence of minerotrophic conditions in our "ombrotrophic" category. Sites characterised as ombrotrophic included raised bogs, blanket bogs, and the extensive bog systems of western Siberia (Kremenetski et al., 2003).

Chronology and age modelling
All sites were 14 C dated using selected aboveground plant remains, except for site 68 where bulk peat was 14 C dated. We recalibrated all the dates from the original studies. For modern (post-AD 1950) 14 C dates, the NH1 postbomb calibration curve was used (Hua and Barbetti, 2004). Remaining dates were calibrated using IntCal09 . Age models for the temporal analysis were based on the program "Bacon", a flexible Bayesian age-depth modelling approach that uses prior information on plausible accumulation rates and their variability and autocorrelation over time Christen, 2005, 2011). Peat cores were divided into contiguous 2 cm segments, and linear accumulation rates were calculated for all individual segments sequentially down the core. Age models were developed based on several million iterations, followed by thinning to remove any autocorrelation between individual model runs, yielding ca. 5000-8000 iterations for each site (Fig. 2a).

Spatial analysis of carbon accumulation
Total carbon accumulation over the 1000 yr was estimated based on a data set of the 90 radiocarbon-dated peat profiles (Tables 1 and 2). The post-1000 yr carbon pool is the difference between total carbon additions from photosynthesis and cumulative respirative carbon release over this interval, reflecting carbon sequestration at a site. We analysed the relationship between total carbon accumulation over the last 1000 yr and climate parameters using a 0.5 • grid, derived from the CLIMATE 2.2 data (Kaplan et al., 2003). Climate parameters included growing degree days above 0 • C (GDD0), cumulative photosynthetically active radiation during the growing season (PAR0), PAR over the growing season, growing season length (days) and the moisture index P /Eq, where P is annual precipitation and Eq is annually integrated equilibrium evapotranspiration calculated from daily net radiation and temperature (Prentice et al., 1993). PAR was calculated from latitude and sunshine hours (Prentice et al., 1993;Harrison et al., 2010).

Temporal variation in carbon accumulation
A composite carbon accumulation curve was constructed based on the subset of 24 well-dated, high-resolution sites with continuous records for the past 1000 yr ( Table 2). The Bayesian age-depth models allowed chronological uncertainty to be included in carbon accumulation curves (Fig. 2b). Table 2 . Characteristics of the high-resolution sites used in the analyses. BDsbulk density sample size; BDibulk density increment depth; dates (other): Pblead 210, Sspheroidal carbonaceous particles, Ttephra, A -Ambrosia pollen rise; scaling for carbon estimates based on C (carbon analyses), LOI (loss on ignition), or 50 % of dry mass.    Christen, 2005, 2011). (b) Carbon accumulation derived from age-depth models, and bulk density and C measurements. Curves are fitted to each of 10 000 possible age models based on Bayesian analysis. Points represent individual samples on different age models, and grey lines are fitted curves for individual models. Only curves fitted incorporating long-term decay and ecosystem maturity (Yu et al., 2003) are shown here.
All age depth models were converted to carbon accumulation using bulk density and carbon or LOI measurements.
We derived different estimates of variability in carbon accumulation rates based on different assumptions about autogenic processes of long-term decay (Clymo, 1984) and ecosystem maturity (Yu et al., 2003). Carbon accumulation rates calculated from our age-depth models and carbon density do not take account of autogenic peat accumulation processes, most importantly the effect of long-term decay. Dead plant material decays rapidly in the surface layers, as the most labile organic matter is broken down quickly by microbial activity. Decomposition rates are much slower (though not zero) in the permanently saturated zone, which contains more recalcitrant organic matter (Clymo, 1984;Belyea and Baird, 2006). If productivity and decay are constant, measured apparent accumulation rates will be higher for more recent peat, and the long-term carbon storage will appear to increase. We accounted for this ecological process by fitting decay curves to each profile (Clymo, 1984). We also tested the effect of "ecosystem maturity", that is the slowing of peat growth under stable conditions because of autogenic limits on the height of the peat surface (Yu et al., 2003). We excluded carbon accumulation changes in the uppermost peat (conservatively approximated here as peat formed after 1850) where relatively rapid aerobic decay is still taking place. We used AD 1850 for this because this is likely outside of the aerobic decay zone for all cores.
The changes in accumulation rates for each site were expressed as differences between observed accumulation and those derived from three models: (1) linear decay model (i.e. no autogenic processes); (2) the Clymo model, which includes long-term decay only (Clymo, 1984): where M is the accumulated carbon, P c is the peat added to the catotelm each year (g C cm −2 ), a c is the catotelm decay constant and t is time; and (3) the extended peat accumulation rate (ExtPAR) model (Yu et al., 2003), which includes long-term decay and ecosystem maturity: where the parameters are the same as those listed above, with the addition of b c , a coefficient that allows the accumulation rate to be modified. Each curve fitting exercise produces estimated values for P c , a c and b c . There are several other more complex models that could be applied to account for longterm decay, but it is often difficult or impossible to determine the most appropriate one, given the subtle variations in the carbon accumulation curves (Belyea and Baird, 2006). Our intention here is to test whether observed variations in the raw carbon accumulation data could be explained by longterm decay and ecosystem maturity. Decay models were fit individually to each carbon accumulation curve derived from the Bacon routine. For the linear model, the fitting was carried out using ordinary leastsquares. For the Clymo and ExtPAR models, optimization was carried out using an iterative orthogonal sampling technique that samples the entire parameter space, then uses a least squares fit to obtain a subset of the parameter space. This subset is then sampled in the next iteration to produce Steps involved in deriving a non-autogenic accumulation curve from a single age profile after fitting of the peat accumulation model with long-term decay and ecosystem maturity (Yu et al., 2003). See text for details. an increasingly well defined parameter space. All following analyses were applied to results from all three versions (Fig. 3).
The number of age models varied by over an order of magnitude between sites (from ca. 2300 to ca. 44 000). To avoid biases towards sites with a higher number of age models, a single age model was randomly sampled from each site. The accumulation rates were interpolated onto a regular time step by taking the median value in a moving window (halfwindow of 25 yr) with a step of 10 yr to avoid bias toward sites with higher sampling resolution. Finally, a time series of accumulation rates was calculated as the median of the 24 interpolated accumulation rates (one per site). This was repeated 10 000 times to provide a fair sampling of the available age models, and gave a matrix of median time series on a regular time step. Finally, this matrix was used to calculate the median and percentile values of the accumulation rates for each time.
The effect of this Monte Carlo resampling of the possible age models (and associated accumulation rate curves) is to give greater weight to the sites with the best-constrained chronologies. In each iteration of the resampling, we took one age model and set of accumulation rates from each site. Sites that are well-constrained will provide age models that are similar in each iteration, and poorly constrained sites will provide age models that are widely different. The end result of this is that well-constrained sites will effectively have a greater weight in the overall composite.
To avoid any bias toward sites with generally very high accumulation rates, a second composite was made, based on transformed values. This followed the methodology used previously for charcoal data (Marlon et al., 2008): (1) minimax transformation of the original accumulation rate time series; (2) Box-Cox transformation to normalise the time series; and (3) z-score calculation. The composite z-scores were estimated using the same procedure as for compilation of the untransformed values described above. The final composite curves are shown in Fig. 4.

Spatial relationships between carbon accumulation and climate
Warming would be expected to increase net primary productivity (NPP) in high-latitude ecosystems because of increased growing season length. The growing season for northern peatlands is appropriately defined as the period of the year with air temperatures above freezing, because bryophytes begin photosynthesis at this threshold, and are the dominant peat-former in most of our sites. PAR, determined by latitude and cloudiness, is the driver of photosynthetic carbon fixation and may also be an important control on NPP. However, higher temperatures could also increase peat decomposition rates through accelerated microbial activity (Ise et al., 2008;Dorrepaal et al., 2009). Linear regression of total carbon accumulated over the last 1000 yr (C) against PAR0 yielded the strongest relationship: with an R 2 of 0.33 (Fig. 5a). In single-predictor regressions, C showed a weaker relationship with GDD0 (R 2 = 0.13, Fig. 5b) and no significant relationship with P /Eq (P = 0.19, Fig. 5c). Residuals from Eq. (3) showed no systematic relation to either GDD0 or P /Eq and inclusion of these additional predictors in a multiple linear regression yielded nonsignificant regression coefficients. The correlation between PAR0 and GDD0 is high (0.83), owing to the growing season length that is shared by both variables. We checked the influence of two apparent outliers with higher PAR0 values on our conclusions. These are the two southernmost sites from Dhakuri (India) and Pinhook (USA). Removing these two sites does not affect the significance of the relationship between PAR0 and 1 ka C (P < 0.0001) but changes the R 2 values from 0.33 to 0.24 and slightly changes the slope from 0.0055 to 0.0049. Thus, it still explains more of the variation than GDD0. The influence of these two sites is not insignificant, but removing them does not impact our main conclusions concerning PAR. Without the two "outliers" total C still shows a positive significant relationship (P < 0.001) with GDD0 but with a change in R 2 from 0.18 to 0.13 and a  Our analyses thus show that total carbon accumulation over the past 1000 yr is linearly related to contemporary PAR integrated over the growing season (PAR0) (Fig. 5a) and that this relationship is stronger than that with growing season warmth expressed as accumulated temperature (GDD0, growing degree days above zero) (Fig. 5b). The implied effect of a warmer climate is to increase NPP to a greater extent than decomposition, suggesting a negative climate feedback in peatlands. Because no relationship was found between total carbon accumulation and moisture, we infer that although an adequate moisture supply is necessary for the presence of peat, above a threshold of moisture availability the effect on carbon accumulation is secondary relative to growing season temperature and light conditions, at least for millennial averages over large spatial scales. Some of the unexplained variability in carbon accumulation probably reflects local hydrological factors, not captured by the large-scale moisture index, as well as other local controls. Peatlands do not occur, for example, in climatically suitable locations with steep slopes. Local topographic and drainage features, as well as internal dynamics, create heterogeneity in peat accumulation that is not represented by our data. Furthermore, our sampling is necessarily biased toward peatlands that exist today, and we therefore cannot establish the threshold for cessation of carbon accumulation resulting from a reduction in moisture balance. Despite these caveats, it is clear that the changes in moisture balance are unlikely to be an important control on peat accumulation during the recent past because of the complete lack of any relationship with macroclimate. If moisture does not play a role in determining peat accumulation rates, this implies that the net balance between NPP and decay is similar under varying hydrological conditions. In dry sites high decay rates must be offset by high NPP, and similarly, wetter sites must have low NPP together with low decay because of the anaerobic conditions.
The results suggest that spatial variability in peatland carbon accumulation over the last 1000 yr is primarily controlled by spatial variation in NPP, which in turn is driven by growing season length (related to temperature) and growing season PAR. There is a statistically significant correlation between large-scale spatial variability in Sphagnum moss productivity and PAR0 that supports this hypothesized mechanism (Loisel et al., 2012). PAR0 incorporates both a PAR and temperature effect through growing season length; we therefore tested the relationship between total carbon and mean PAR over the unfrozen season. The reduced R 2 (0.13, p < 0.0001) compared with PAR0 (R 2 = 0.33), suggests that PAR has an effect independent of temperature. A much weaker relationship with growing season length (R 2 = 0.08, p < 0.01) implies growing season length is of subsidiary importance. Taken together, these results support our hypothesis that peatland carbon accumulation is driven by PAR over the growing season.

Temporal changes in carbon accumulation
There is an overall downward trend in the composite carbon accumulation rates from AD 1000 to 1850 (slope −0.0026 g C m −2 yr −2 , p < 0.0001), implying reduced peat accumulation during the LIA (Fig. 6a). The decline appears to be greater in the latter (post AD 1400) period than the earlier part of the record. The downward trend is present in the raw data as well as in the curves that include long-term decay and ecosystem maturity at individual sites, showing that the direction of change is insensitive to any assumptions concerning long-term ecosystem processes (Fig. 4). The magnitude of this change differs depending on the time period chosen. Here we compare long-term median values for the periods 1000-1425 and 1425-1850 AD, broadly corresponding to the times used to define the MCA and LIA (Mann et al., 2008(Mann et al., , 2009Jansen et al., 2007). We also compare the difference at the start (1000) and end (1850) of this period based on the slope of a regression line through the data. Using both approaches means that the analyses are robust to short-term fluctuations in the data. If long-term decay is taken into account, the difference between the mean of the median accumulation rates between these two periods is 2.43 (± 0.92) g C m −2 yr −1 . The estimate for this is smaller (1.01 ± 0.89 g C m −2 yr −1 ) if some of the decay is compensated for by ecosystem maturity at individual sites (Fig. 4). The difference between the start (1000) and end of the data (1850) is 6.05 ± 5.40 and 3.15 ± 5.08 g C m −2 yr −1 , respectively, for a decay only and decay plus ecosystem maturity model (Fig. 4). Shorter term changes in the carbon accumulation curve may be related to temperature or other climate variables, but the changes are not sufficiently robust to draw firm conclusions about changes in carbon accumulation on (sub)centennial timescales.

Climate controls on carbon accumulation
One explanation for the reduction in carbon accumulation between the MCA and LIA is that decreased temperatures reduced NPP through shorter growing seasons (reduced GDD0). Shorter growing seasons would also reduce accumulated PAR0 assuming light levels remained unchanged. We tested this hypothesis by using the spatial relationship between total 1000-yr carbon accumulation and climate variables shown in Fig. 5 to calculate the effect of MCA-LIA cooling from palaeoclimate records (Fig. 6b) on carbon accumulation and comparing it with the observed changes in carbon accumulation (Fig. 6a).
We used the IGBP-DIS soil carbon gridded data set (http://daac.ornl.gov/SOILS/guides/igbp-surfaces.html) to select all grid cells occupied by northern peatlands, and summed their carbon accumulation rates as predicted by PAR0 from Eq. (1) (multiplied by the grid cell areas) to estimate the total carbon sink in northern peatlands. "Peatlands" were defined as 0.5 × 0.5 • grids north of 40 • N that contain 10-min IGBP soil C grids that are all > 31 kg C m −2 (Wania et al., 2009). This is conservative, designed to focus on the biggest peatland areas that dominate the global peatland C cycle. PeatStash (Gallego-Sala et al., 2010) was used to calculate the accumulated PAR0 by summing the daily PAR0 over the growing season (days above freezing) for each peat-land grid cell. The daily PAR0 is obtained by integrating the instantaneous PAR between sunrise and sunset (Harrison et al., 2010). The seasonal accumulated PAR0 depends on latitude and cloudiness, and indirectly on temperature, because the temperature determines the length of the growing season, i.e. which days are included in the seasonal accumulated PAR0 calculation. We calculated changes in GDD0 and PAR0 that would result from the change in temperature inferred from the palaeoclimate reconstructions over the last 1000 yr (Fig. 6b). The temperature difference in the palaeoclimate records was based on calculations similar to those made for carbon accumulation over the same time periods. The median temperature difference between the two periods 1000-1425 and 1425-1850 is 0.116 ± 0.02 • C. We calculated the influence of this change in temperature on the duration of the growing season (Fig. 7a), and applied this to the relationships between carbon accumulation, GDD0 and PAR0 derived from the spatial analysis (Fig. 5a, b). Changes in PAR0 and GDD0 for each peatland grid cell were calculated by adding 0.116 • C to the input climate data set used in PeatStash.
The results suggest that the sensitivity of GDD0 and PAR0 to temperature is too small for either of these to provide the sole explanation for the observed change in carbon accumulation rate over the last 1000 yr. Over the peatland areas, a 0.116 • C increase results in a mean change of +18.5 GDD0 (5-95 % range, 12.3-24.7 GDD0), which is predicted by the regression equation to result in a mean increase of 0.16 g C m −2 yr −1 (from 0.11 to 0.21 g C m −2 yr −1 ). PAR0 increases by 13.5 mol photons m −2 season −1 (range 0-37.4 mol photons m −2 season −1 ), which is predicted by the regression equation to result in an increase of 0.07 g C m −2 yr −1 (from 0 to 0.21 g C m −2 yr −1 ). If the magnitude of temperature change estimated from the Northern Hemisphere records is assumed to be correct, it is therefore unlikely that the observed changes in carbon accumulation of 1.0-2.4 g C m −2 yr −1 are a result of temperature and growing season changes alone. The conclusion is similar if the magnitude of temperature change over the period 1000-1850 is used for this calculation. In this case, a temperature change of −0.266 • C is associated with a reduction in carbon accumulation of 3.2-6.1 g C m −2 yr −1 . However, mean changes in GDD0 would be 43 GDD, and 31 mol photons m −2 season −1 in PAR0, resulting in predicted carbon accumulation reductions of only 0.37 and 0.17 g C m −2 yr −1 , respectively. It is possible that Northern Hemisphere averages underestimate temperature changes in peatland regions, as estimates of maximum MCA-LIA difference range from > 0.3 to as much as 1.8 • C in a few locations (Mann et al., 2009). Repeating the analysis with a larger 1 • C change suggests a response of 0.4 to 0.9 g C m −2 yr −1 from PAR0, which is closer to but still less than the 1.0-1.4 g C m −2 yr −1 reduction in carbon accumulation shown in the data over the MCA-LIA.
If temperature change alone is not the driver of carbon accumulation changes, an additional hypothesis is that reduced PAR from an increase in growing-season cloudiness significantly reduced NPP during the LIA, causing a fall in peatland carbon accumulation. We therefore also examined the sensitivity of PAR0 and carbon accumulation to changes in cloudiness (Fig. 7) by altering, in addition to 0.116 • C warming, the annual sunshine hours by a range of +/− a percent of the input climate data set. Sensitivity tests (Fig. 7a) show that a 4 % increase in sunshine hours on top of the influence of a +0.116 • C change (equivalent to 199 mol photons m −2 season −1 ) could result in an average change of +1.1 g C m −2 yr −1 over the peatland areas. These are averages; the change in PAR0 has a 5-95 % range of 137-255 mol photons m −2 season −1 and change in carbon accumulation between 0.8 and 1.4 g C m −2 yr −1 for this sensitivity analysis. A geographical pattern emerges in these simulations where greater sensitivity of PAR0 to temperature and cloudiness occurs in lower latitudes (Fig. 7b). The existence of a positive relationship between PAR0 and peat carbon accumulation is further supported by Fig. 7a, which shows that high soil carbon grid cells in the independent soil carbon density data are located in areas with higher PAR0. Our finding that spatial patterns of carbon accumulation can be explained by spatial variability in climate, specifically PAR0, implies that the temporal variations in carbon accumulation over the last millennium may also be explained by the same climatic variables. The data are consistent with this space for time substitution argument in that the warmer MCA has higher rates of carbon accumulation than the cooler LIA, i.e. longer, warmer growing seasons in the MCA increased carbon accumulation in comparison to the LIA. That the MCA was warmer than the LIA is not contested (Jansen et al., 2007), but we also suggest from our sensitivity analysis that temperature alone cannot explain the magnitude of the change in observed carbon accumulation over the last millennium. Changes in PAR received by the plants due to cloudiness or some other influence such as snow cover (increased depth/snow lie after T > 0 during the LIA) or diffusivity are also required to generate the observed changes, assuming that temporal sensitivity is the same as that derived from spatial relationships. We explored the potential reduction of PAR by late-lying snow cover by calculating the PAR for days without snow and with temperature > 0 • C. The correlation between snow-free PAR0 and total carbon accumulated over the last 1000 yr was not significant (p = 0.23), so cloudiness is the more likely cause of change in PAR and carbon accumulation. The idea that the Northern Hemisphere LIA was characterised by greater summer cloudiness is consistent with historical documentary data (Grove, 2004), although recent findings from δ 13 C in Fennoscandian tree ring studies suggest that there may have been regional differentiation of this tendency (Young et al., 2010;Gagen et al, 2011).

Changes in the strength of the peatland carbon sink
An average decrease of 1.01 to 2.43 g C m −2 yr −1 in carbon sequestration rates in peatlands between the MCA and the LIA represents a reduction in the strength of the sink of 0.0035 to 0.0085 Gt C yr −1 , assuming a total northern peatland area of 350 million ha (Gorham, 1991;Tarnocai et al., 2009). This should be compared with our estimate of the median net accumulation rate between 1000 and 1850 of 26.1 (sd. 1.47) g C m −2 yr −1 or 0.091 (sd. 0.005) Gt C yr −1 . This latter figure is similar to the average figure of 0.096 Gt C yr −1 (29 g C m −2 yr −1 ) estimated for the whole of the Holocene (Gorham, 1991). Neither of these figures for long-term net accumulation take into account the very slow long-term decay of the deeper peat layers. Decay of this store might reduce the actual average net ecosystem carbon balance to around 0.042 Gt C yr −1 (Yu, 2011) for the last millennium, depending on how much peat had accumulated before AD 1000, its permafrost status, and its average age and long-term decay rate. Thus, a reduction in carbon sink of 0.0035 and 0.0085 Gt C yr −1 represents a decrease of about 8 % and 20 % in net ecosystem carbon balance of northern peatlands. This is an average figure for the two 425-yr periods chosen to represent the MCA and LIA here. The extremes in change are much greater; for example, the difference between the 25-yr mean at 1250 and 1850 is 5.9 ± 5.2 to 3.9 ± 5.1 g C m −2 yr −1 , or a 33 to 45 % reduction in the strength of the CO 2 sink.
The main decrease in atmospheric CO 2 concentrations of approximately 6 ppm occurred over the period 1150 to 1750 (using the average curve in Fig. 6c). For this time period the decline in peatland carbon sink is between 2.79 (± 4.88) and 0.77 (± 4.87) g C m −2 yr −1 , equivalent to 0.0027 to 0.0097 Gt C yr −1 if applied to all northern peatlands. Assuming a linear decline in strength of the peatland sink over this period, this would represent a change equivalent to about 0.38 (± 2.41) to 1.38 (± 2.41) ppm in the atmosphere. Although this is a small change, because the observed decline in atmospheric CO 2 is unequivocally attributed to increased land carbon storage by the concomitant increase (Elsig et al., 2009) in the δ 13 C of atmospheric CO 2 , the reduced carbon sequestration in northern peatlands has to be added to the amount of CO 2 reduction to be explained by carbon uptake, presumably as a result of suppressed decomposition of soil organic matter or increased NPP, in other terrestrial ecosystems.

Conclusions and implications
We have shown that there has been a small negative feedback to climate from changes in Northern Hemisphere peat accumulation over the last 1000 yr. The direction of the peat accumulation-climate feedback is supported by large-scale spatial patterns of peat accumulation over the past 1000 yr in relation to modern climate gradients. Northern peatlands sequestered carbon at a higher rate during the MCA than during the LIA. Although the magnitude of change is small (approximately 1 ppm CO 2 over the MCA-LIA transition), this suggests that carbon accumulation in northern peatlands may also change in response to future anthropogenic climate warming. Our spatial analyses indicate that the strength of any feedback effect could depend on changes in cloud cover as well as changes in temperature and growing season length. Continued carbon accumulation depends on sustained adequate moisture availability to maintain peatland growth. The sites we sampled all have adequate moisture supply at present and probably also during most of the last 1000 yr, but this may not continue in the future. Climate model projections suggest that most of the high latitudes will experience higher summer temperature and higher precipitation, with decreases in soil moisture for some regions (Meehl et al., 2007). Midlatitude peatlands in locations such as western Europe are most vulnerable, especially as summer precipitation is projected to decrease. If the threshold between presence and absence of peatlands on the gradient of precipitation/equilibrium evapotranspiration ratio is crossed (Fig. 5c), some peatland areas (including some of the blanket bogs characteristic of extremely oceanic climates; Gallego-Sala et al., 2010;Gallego-Sala and Prentice, 2013) may stop growing. Changes in peatland extent in response to climate change occurred in the past, including during the last millennium (Finkelstein and Cowling, 2011), and will occur in the future (Gallego-Sala and Prentice, 2013). Our analysis suggests that a reduction in peatland area in the midlatitudes could be compensated by increased carbon accumulation across the very large areas of peatland in higher latitude regions, but only if the majority of these peatlands retain sufficient moisture and there is no significant increase in cloud cover. Current projections of cloud cover are highly uncertain for the northern peatland regions; some areas are projected to have less cloud and others more cloud (Meehl et al., 2007). A further compensating factor for peatland loss is that peatlands may increase in extent in high-latitude areas that are currently too cold and dry for peat formation. Other major uncertainties concerning feedbacks between peatlands and climate change still exist, including changes in fire regimes, nitrogen deposition, permafrost thaw and the role of methane emissions. Our analyses show that only about a third of the variability in C accumulation is explained by PAR and other factors must play a role. We might expect climate changes at some of the sites used in our analysis to move beyond the climate envelope explored here, and this could affect their contribution to the carbon budget. However, based on our analyses of carbon accumulation over the past millennium, and contrary to the conclusions from soil decay models (Ise et al., 2008;Dorrepaal et al., 2009), we suggest that carbon sequestration may increase in many high-latitude peatlands in response to future climate warming over the next century.