Greenhouse gas balance of a semi-natural peatbog in northern Scotland

Northern peatlands have been accumulating organic matter since the start of the Holocene, and are now a substantial store of terrestrial carbon. However, their current status as carbon sinks is less clear, because of the possible effects of climate change, air pollution, grazing and drainage etc., and the difficulties of accurate measurement with suitable time resolution. Such measurements are particularly lacking in the UK. Here, we present multi-year eddy covariance measurements of the carbon fluxes at a relatively undisturbed ombrotrophic blanket bog in the Flow Country of northern Scotland. The site consistently acted as a moderate sink for CO2 over all the measurement years (mean net ecosystem exchange (NEE) of −114 g C m−2 y−1), similar in magnitude to other measurements in the boreal and tundra zones, and rather higher than the existing measurements at other sites in the UK and Ireland. Generally, the NEE of CO2 was relatively insensitive to moderate inter-annual variations in weather. Non-CO2 losses comprised 11% of gross primary production, mainly from methane emissions. Accounting for these terms, the net ecosystem carbon balance was −50 g C-CO2 eq m−2 y−1. The contemporary carbon sink was larger than estimates from local peat cores, based on peat accumulation over the last several thousand years, but in the middle of the range of estimates which used spheroidal carbonaceous particles to estimate peat accumulation rates over the last century.


Introduction
In Scotland, peatlands are estimated to hold 1.62 Pg C, 56% of the total soil carbon (Chapman et al 2009), and have the potential to act as a major source or sink for carbon. Globally, northern hemisphere peatlands have been accumulating carbon since the last glacial maximum, at a rate of around 19 g C m −2 y −1 (Yu 2011), resulting in a contemporary stock of around 547 Pg C. This carbon sink has varied with long-term variations in climate, and is thought to have declined slightly from the Medieval Warm Period to the Little Ice Age (ca. AD 1550-1850, Charman et al 2013). As well as climatic change, peatlands in the UK have been subject to elevated CO 2 , air pollution and land management practices including grazing, burning, drainage and afforestation. However, because of the lack of longterm research and the difficulties of estimating carbon uptake, the current status of peatlands in Scotland as carbon sinks is unclear.
Because of the limits on depth/time resolution, current trends are difficult to discern from peat core dating methods. Measuring change in carbon stocks directly, by sequential measurements of peat depth, bulk density and carbon content in time, is generally not possible because of the huge sampling difficulties and the time periods involved (Chapman et al 2013). As a means of estimating the current carbon sink, we present here direct CO 2 flux measurements using a micrometeorological approach (eddy covariance, Baldocchi 2003) over a six-year period at a near-pristine peatbog in northern Scotland. Additional measurements of methane CH 4 emissions were made to establish the net greenhouse gas (GHG) balance, and including estimates of fluvial losses from the catchment.
The site is in the 'Flow Country' of Caithness and Sutherland, an extensive area of Sphagnum-dominated blanket bog, characterised by pool systems, with a mild, wet, oceanic climate (Lindsay et al 1988).
Blanket bog in Britain develops as a response to cool, oceanic conditions, distinct from tundra associated with permafrost environments and boreal regions which harbour most of northern peatlands. In the 'Flow Country' of northern Scotland, the topography, geology, and oceanic climate combine to produce an environment where blanket bog development reaches an extreme, with a floristic composition which is globally unique. The area has the largest expanse of deep peats in the country (Chapman et al 2009), and constitutes a large proportion of the total peatland area (and carbon stock).
Although there are several comprehensive studies in North America and Fenno-Scandinavia (Roulet et al 2007, Nilsson et al 2008, few direct measurements are available with which we can quantify the contemporary GHG balance of UK peatlands. There are, in fact, few (if any) published long-term ecosystem-scale measurements of CO 2 fluxes over blanket bog in the UK (Billett et al 2010). Short-term micrometeorological measurements of CO 2 and CH 4 fluxes have been made in the Flow Country (Beverland et al 1996, Hargreaves andFowler 1998), but the duration of measurements was limited to a few weeks. Long-term measurements have been maintained for over ten years at Auchencorth Moss in central Scotland (Helfter et al 2014), but the tower site represents a mix of acid and marshy grassland and modified bog (UK National Vegetation Classification (NVC) types U2, U4, U5, MG10 and M15, Rodwell 1998), rather than a typical semi-natural blanket bog (NVC M17, M18, M19). Virtually all other data on CO 2 and CH 4 fluxes from UK peatlands come from chamber methods. Because they measure over only small areas (typically 0.1 m −2 ) and on relatively few occasions, there are real limitations in upscaling to the whole ecosystem over annual time scales. Despite the current interest in the possible benefits of peatland restoration as a means of sequestering carbon, and thereby mitigating climate change, we still lack the basic data on contemporary peatland carbon sequestration in the absence of human intervention. The measurements presented here provide the first long-term data on GHG balance at a near-pristine peatland in the UK, as a baseline against which restoration activities might be compared. The aims were: to establish the GHG balance over multiple years, including estimates of all the major budget components; to relate the measured fluxes to driving variables, and; to compare the measured carbon balance with estimates derived from peat cores in the local area.
2. Field site and methods 2.1. Field site Micrometeorological equipment was installed in April 2008 at a site on the Royal Society For Protection Of Birds nature reserve at Forsinard in Sutherland, Scotland. The reserve covers an area of 215 km 2 and ranges in elevation from 44 to 580 m above sea level, with most of the deep peatlands between 120 and 438 m elevation. The tower site is to the south-west of the Cross Lochs 58 22 13 N, 3 57 52 W ( )  ¢   ¢  , on a large, flat expanse of blanket bog with pool systems typical of Caithness and Sutherland Flows. There was a fetch of at least 3 km over blanket bog in the direction of the prevailing wind (south-west), with larger extents to the south and east. Approximately 1.5 km to the west was an area of plantation forest, while 600 m to the north was an area of plantation forest which had been felled five years previously. The sensors were installed at a height of 3 m, such that the footprint of the measurements would extend over an area of bog with minimal influence from any forest areas. Under turbulent conditions, the source area contributing to the measured flux was generally within a radius of 100 m from the tower. The bog area itself had not undergone any management practices, such as burning or draining, and was free from any direct human intervention as far as we could tell. Air pollution is relatively low in Northern Scotland, with current N deposition rates of less than 5 kg N per hectare per year (Fowler et al 2005). Grazing by red deer (Cervus elaphus L.) may have some influence on the bog, but deer numbers on the open bog were low. Given the vegetation composition, the high water table, and intact pool systems characteristic of natural bog in the area, the site is potentially as close to a pristine state as can be found in the UK.

Carbon dioxide exchange
A micrometeorological approach, eddy covariance, was used to make near-continuous measurements of the surface fluxes of carbon dioxide (CO 2 ) and water vapour (Baldocchi 2003). In brief, the net ecosystem exchange (NEE) of CO 2 , F c , is given by: where w′ is the instantaneous deviation of the vertical windspeed from the mean, and c′ is the instantaneous deviation of the CO 2 concentration from the mean. In order to equate the eddy covariance, w c , ¢ ¢ with the surface flux, we need to make assumptions about stationarity, and horizontal homogeneity, and apply coordinate rotation and Reynolds averaging (Lee et al 2006). The three components of windspeed were measured by an ultrasonic anemometer (USA1, Metek GmbH, Elmshorn, Germany), mounted at a height of 3 m. CO 2 and H 2 O concentrations were measured by an open-path infra-red gas analyser (IRGA)(Li-7500, Licor Corp., Nebraska, USA) with a response time 40 Hz. A data logger (CR3000, Campbell Scientific, Loughborough, UK) logged the data from these instruments at 10 Hz. The Li-7500 IRGA was mounted at a zenith angle of 15°, pointing due north, so as to minimize direct insolation.
Ancillary measurements were made with sensors for solar radiation (Kipp & Zonen, Delft, Netherlands), quantum flux density (Q, in μmol m −2 s −1 ) (SKP215, Skye Instruments, Llandrindod Wells, UK), air temperature (T air in°C) and relative humidity (Model CS215, Campbell Scientific, Loughborough, UK), soil temperature at 10 cm depth (T soil in°C) (Model 109 T, Campbell Scientific, Loughborough, UK), soil moisture (CS616, Campbell Scientific, Loughborough, UK), water table height (z water , in cm) (PDCR1830, Druck, Burlington, VT., USA) and rainfall (ARG100, Campbell Scientific, Loughborough, UK). Based on T air data for each month, season and year, we calculated thermal time (t thermal in°C · days) as the sum of daily mean air temperature above a 5°C threshold, and the number of potential growing season days (n gt5 , where mean air temperature was >5°C for at least seven contiguous days. Remotely-sensed normalized difference vegetation index (NDVI) data were obtained for the 250 m pixel centred on the tower from the MODIS instrument (ORNL DAAC 2014).
Power was supplied by a wind turbine (Model 910-3 Furlmatic, Rutland, Notts., UK), solar panels, and in later years, a fuel cell. These charged a large array of deep-cycle sealed lead-acid batteries. The datalogger controlled power consumption by switching off the sonic anemometer and the IRGA when battery voltage was low, but maintained the ancillary measurements as far as possible. Interruptions to the system operation were most common in the winter when solar energy was at a minimum, and access to the site was most difficult. Half-hourly data were sent by telemetry via the mobile telephone network to the base institute.
The raw 10 Hz data were processed using the EddyPro ® software, with further analysis in R (R Core Team 2014). In the processing, we applied double coordinate rotation (vertical and crosswind), spike removal, block averaging, and time lag removal by covariance maximization. Correction for the frequency response of the system, both high and low-frequency losses, were made using the method of Moncrieff et al (1997). The Webb-Pearman-Leuning correction for density fluctuations was applied on a half-hourly basis (Webb et al 1980). The correction for an additional sensible heat flux in the path of the IRGA (Burba et al 2008) was not applied because of the off-vertical mounting, the generally low insolation and mild temperature environment (but see section 4).
Gaps in the flux data were filled using the algorithm of Reichstein et al (2005), as implemented in the R package REddyProc (Reichstein and Moffat 2014). Longer gaps were filled by fitting a general additive model to the data, using the R package mgcv (Wood 2006). Partitioning of NEE fluxes into gross primary production (GPP) and ecosystem respiration (R eco ) was also calculated using the algorithm of Reichstein and Moffat (2014). The uncertainty in observed half-hourly fluxes was estimated using the method of Finkelstein and Sims (2001). For gaps in the flux data, uncertainty in the predicted fluxes was calculated from the algorithm of Reichstein et al (2005) and from the prediction intervals of the general additive model.

Methane exchange
Methane and nitrous oxide fluxes were measured by the static chamber method (Hutchinson and Mosier 1981). Cylindrical PVC collars (either 63 cm in diameter and 29 cm high, or 38 cm in diameter and 25 cm high in later years) were inserted into the soil and left in place for a number of months or years. Twenty chambers were located in the immediate vicinity of the flux tower, mainly in the north-west sector, so as not to disturb the primary fetch. Further chambers were located on other relatively undisturbed bog sites in the surrounding area, so as to cover a range of spatial variabilty, although these measurements preceded the timespan of the eddy covariance measurements (2003)(2004)(2005).
On each sampling occasion, a lid was sealed on top, and left in place for 30-40 min. A 5 V fan was run at a low rate inside the chamber, to ensure sufficient mixing. Four 20 ml samples were removed by syringe through a 3-way tap or rubber septum, stored in vials or tedlar bags, and analysed on a gas chromatograph (5890 series II, Hewlett Packard), together with replicates of three or four standard gases with known concentrations. For each sequence of gas samples from a chamber, the flux was calculated as: where F is gas flux from the soil (μmol m −2 s −1 ), C t d d 0 is the initial rate of change in concentration with time in μmol mol −1 s −1 , ρ is the density of air in mol m −3 , V is the volume of the chamber in m 3 and A is the ground area enclosed by the chamber in m 2 .
The parameter C t d d 0 was calculated using linear and non-linear asymptotic regression methods (Levy et al 2011). Using a mixture of goodness-of-fit statistics and visual inspection, the regression method that provided the best fit for the time series of concentration was chosen for each individual measurement. With this method of flux calculation, any non-linearity should be accounted for as far as possible. However, the time resolution (approximately 10 min) limits the detectable degree of non-linearity in the initial concentration change, so there remains some potential for underestimation of fluxes (Cowan et al 2014). Fluxes of nitrous oxide were consistently below the limit of detection of the measurement system and only methane fluxes are presented here. Figure 1 shows the time series of the main environmental variables over the six-year period. The site has a mild oceanic climate, with daily mean T air never exceeding 20°C in summer, and few extended periods below 0°C in winter. This contrasts with most peatland areas in the northern hemisphere, which have extended freezing periods and warmer summers. Rainfall was spread relatively evenly across the six years, though with notable dry spells in the summers of 2008 and 2010. These show clearly in the record of water table height. Quantum flux density was very consistent from year to year. NDVI shows a strong seasonal signal, lagging behind that of Q because of the delay introduced by plant growth. Some of the departure from the seasonal pattern will be instrumental noise, but there are plausible real signals in NDVI, such as the drier periods in 2008 and 2010, and the warm spring of 2012. Figure 2 shows the observed NEE over the same period, and the gross fluxes calculated from this. The data show the seasonal cycle clearly, but clear differences between years are not very apparent, especially given the variabilty in the data and gaps in the time series. Detailed analysis of some periods with more extreme conditions (e.g. summer 2008) did not reveal clearer relationships between NEE and water table height or wind direction. The highest NEE values occurred in summer 2009, when the environmental variables were not particularly different from the average.

Statistical analysis
We analysed statistically the importance of different driving variables over a range of time scales, from daily to seasonal, using multiple linear regression. With data of this nature, it is difficult to separate environmental effects because of their covariation (e.g. radiation and temperature are strongly correlated). However, models are still useful for summarising the results and making empirical predictions within reasonable limits.
One approach is to start with a wide set of variables and use stepwise regression to identify 'optimal' models, according to the Aikake Information Criterion (AIC). We note the limitations of this approach (e.g. see Whittingham et al (2006)), and that this provides only a crude indication of the more important factors. To keep models simpler, we generally excluded interaction terms, with the exception of the Q × NDVI interaction term, which we included in the initial model for each stepwise procedure, because of the theoretical expectation that these will have a multiplicative effect. On a daily time scale, the optimal model for NEE included the terms Q, NDVI, Q × NDVI, T air , T soil and z water , and explained 28% of the variance in NEE (quantified by the adjusted r 2 throughout). The residual standard error (RSE) was 1.02 on 457 degrees of freedom (df). Statistical analysis was applied only to observed NEE data (i.e. without gap-filling) which met the highest quality control conditions, and where there was at least 75% coverage in any averaging period.
On a monthly time scale, a strong relationship between NEE and Q is clear (figure 3), along with lesser relationships with measures of temperature and NDVI. The optimal model comprised the terms Q, NDVI, Q × NDVI, T air , thermal time (t thermal ) and z water , and explained 83% of the variance in NEE (RSE: 0.275, 19 df).
We also analysed the data averaged over the four seasons as conventionally defined in meteorology (winter = DJF, spring = MAM etc). We only analyse GPP and R eco at seasonal scale, so as to average out short-term random error which arises in the flux partitioning process. At seasonal scale, we also included the effects of air temperature and vegetation state (as measured by NDVI) in the preceding two seasons (denoted T t − 1 and T , and NDVI. Again, this explained 92% of the variance in GPP (RSE: 0.281, 13 df). Variability in R eco was harder to explain; the best model had a lower adjusted r 2 , of 0.50, and included Q, NDVI, Q × NDVI, T air , rainfall, z water , T soil , T t , and n gt5 (RSE: 0.348, 10 df). On a yearly time scale, there were too few points to fit multiple regression models, and relationships in the data are not obvious ( figure 4).
Another approach to identifying the important explanatory variables is to examine their correlations with the flux variables. Because the flux partitioning algorithm uses Q and T air , to avoid circularity, we remove the effect of these variables, and analyse the remaining variation. To do this, we first fit a multiple linear regression model with Q and T air to each of the seasonally-averaged NEE, GPP and R eco . We then analyse the correlations between the residual variation in NEE, GPP and R eco and the other environmental variables (figure 5). These residuals show reasonably strong correlations with VPD and rainfall, and weak correlations with NDVI and z water . For the other variables, strong correlations are not apparent. The pattern is similar across all three flux variables, except that the sign is reversed for NEE, because of the sign convention used. Compared with the step-wise approach, this identifies a similar but simpler model, with the terms VPD, rainfall, and NDVI, as well as Q, T air specified from first principles. Effects of temperature or vegetation state in preceding seasons are not apparent here. The results were similar when the analysis was applied at a monthly time scale. Partial correlations (calculated following Kim (2012)) were also examined, and can in principle remove some of the co-linearity between explanatory variables, but do not show a clearer pattern here. This is probably because of the number of variables and the extent of their covariation. Figure 6 shows the cumulative CO 2 exchange over the measurement period. There is a clear seasonal cycle, but a consistent net uptake in the long term, which averages 114 g C m −2 y −1 (table 1). Random uncertainty is estimated at 9.5 g C m −2 y −1 (95% confidence interval) by accumulating the values produced by the algorithm of Finkelstein and Sims (2001) for each half-hour when observations were available, and from the algorithm of Reichstein et al (2005) when gap-filled. The results suggest that NEE tends to be rather similar across years, despite some fluctuations in climate. Partly, this is because the site is always wet, and the range in water table height is relatively small. Much of the variation in NEE (74%) was accounted for by Q (the basic input to photosynthesis) and T air , and the contribution of less direct factors (such as VPD, NDVI, thermal time etc) is less clear. Figure 7 shows the response of methane flux to water table height and soil temperature at the flux tower site and similar sites nearby. The flux tower site was similar to the other sites, in terms of water table height and methane fluxes. Around 35% of the variance was explained by water table height and temperature, which is not unusual for data of this kind. A substantial fraction was explained by vegetation species composition within the chamber (Levy et al 2012), but this does not provide a useful predictive variable in this context. A simple model was fitted to the data, with a two-order polynomial term for water table height and a linear term for soil temperature, and used to predict methane flux over the whole measurement period. Cumulative values were calculated over the six-year measurement period, and gave a mean annual methane flux of 4.3 g C-CH 4 m −2 y −1 . To express this in units of CO 2 equivalents, we convert to a mass of methane and multiply by 34 (the global warming potential of methane over a 100-year time period including carbon cycle feedbacks, IPCC 2013). From this, we obtain an estimate of 53.5 g C-CO 2 eq m −2 y −1 for the mean annual methane emission (table 1).

Methane exchange and fluvial carbon losses
To calculate the net ecosystem carbon budget, we calculate the gross input of carbon via GPP and losses via ecosystem respiration (using the flux partitioning  algorithm), and we subtract losses from methane emission and fluvial carbon in streamflow. Based on observations of discharge and carbon concentrations in the River Halladale 10 km downstream of our site, Hope et al (1997) calculate fluvial carbon losses from the catchment to be 10 g C m −2 y −1 . A more recent study obtained the same value for the immediate subcatchment of our site (Dinsmore, unpublished data), using the same methodology as described in Dinsmore et al (2013), and which included dissolved organic carbon, dissolved inorganic carbon, aqueous CO 2 and particulate organic carbon components. Putting all the terms together, the site still acts as a net sink for CO 2 of −99.4 g C m −2 y −1 , although the net sink in terms of CO 2 equivalents is reduced to −50.2 g C-CO 2 -eq m −2 y −1 when methane emissions are expressed in terms of their the global warming potential (table 1). The non-CO 2 losses constitute 11.1% of GPP on the latter basis, so are small compared to respiration, but not negligible. Hargreaves and Fowler (1998) previously made CO 2 flux measurements at a site 30 km to the east at Loch More. Although they did not attempt to scale up to an annual CO 2 flux because of the short measurement campaign duration, the pattern and magnitude in their CO 2 flux data is similar to ours, reaching an asymptote of around −6 μmol m −2 s −1 in full light. The annual NEE at our site is very close to the mean for seven northern peatland and wet tundra sites of   Table 1. Greenhouse gas budget terms for the site averaged over the six-year measurement period, in units of carbon (g C m −2 y −1 ) or CO 2 equivalents (g C-CO 2 -eq m −2 y −1 ), accounting for the higher global warming potential of methane in the latter. Net ecosystem exchange (NEE) is the balance of gross primary production (GPP) and ecosystem respiration (R eco ). The net ecosystem carbon balance (NECB) is the balance of NEE minus fluvial losses and methane. The closest site geographically with equivalent data shows a rather lower NEE value (−64.1 g C m −2 y −1 , Helfter et al 2014), but the management history and vegetation at this site is somewhat different, being a mix of acid and marshy grassland and modified bog. More similar is the blanket bog site at Glencar in Ireland Kiely 2010, McVeigh et al 2014), where a 10-year record of NEE data averages −56 g C m −2 y −1 , again lower than our site. To what extent our site is representative of intact bogs in the wider region is difficult to assess without replication. Particularly given the interaction between the spatially patchy distribution of pool systems, which may act as 'hot spots', and the spatio-temporal dynamics of the flux footprint, there will be some 'sensor-location bias' in our data (sensu Schmid and Lloyd 1999).

NEE of CO 2
In one of the longest data sets, Peichl et al (2014) reported a 12-year record of NEE measurements at a boreal fen in Sweden. Similar to our results, they found a general lack of sensitivity of NEE to moderate inter-annual climate variations. Our site was more consistently wet, and we generally do not see negative anomalies of the water table having a strong influence on fluxes, as they found. Our driest year (2008) did have the least negative NEE, but this seems not to be associated with low water table in the summer-time. Indeed, the trend is opposite in our data, with low water table associated with more negative NEE (figure 3), though this may be a spurious correlation with other seasonal effects. We do see an apparent effect of rainfall and VPD on NEE, with higher rainfall and humidity associated with a reduction in CO 2 uptake. This is similar to that reported by Nijp et al (2015), at the Swedish fen site. They suggest the phenomenon is explained by rain-associated reduction in Q, rather than by a direct effect of rain on plant water status. However, here we see the effect after the effects of Q and T air have been accounted for ( figure 5). This suggests we may be seeing an influence on plant physiology as well as the purely climatic correlations between Q, cloud cover and rainfall.
Both Peichl et al (2014) and Helfter et al (2014) report a correlation between pre-growing season T air and summer-time GPP, which presumably arises as winter and spring weather conditions affect the subsequent state of vegetation. Although preceding season T air is chosen as a variable in some of the step-wise model selections, there is no strong evidence for this in the correlation analysis here, once the effects of current Q and T air are accounted for (figure 5). A longer time series would be needed to draw stronger conclusions.
Random uncertainty in the measurements themselves is rather small, and the uncertainty added by filtering outliers and filling gaps in the data with model predictions is considerably larger. However, there is potentially a much larger systematic uncertainty associated with the instrument sensible heat correction (Burba et al 2008, Reverter et al 2011. For example, Oechel et al (2014) found that applying this correction at an Alaskan Arctic tussock tundra site turned an apparent sink of −74 g C m −2 y −1 to a source of +14 g C m −2 y −1 . Even at temperate sites, the correction could potentially lie in the range 140-190 g C m −2 y −1 (Reverter et al 2011), andLund et al (2015) show that this dominates the uncertainty at their Norwegian site. If applied to our data, the correction would give a reduction of 146 g C m −2 y −1 , yielding a NEE of +32 g C m −2 y −1 . Given the offvertical mounting, the generally low insolation and mild temperature environment, and partly for clarity in what we have calculated, we report the uncorrected fluxes here. These are comparable with the bulk of the literature, where the correction has not been applied. To apply the correction properly, one would have to decide on an uncertain temperature threshold below which to apply the correction (Lund et al 2015), and account for the off-vertical geometry of the instrument. The recent development of low-power instruments which do not have the same issue will reveal whether (or to what extent) the correction should be applied retrospectively.

Methane exchange
Our estimate of the annual methane flux of 4.33 g C-CH 4 m −2 y −1 is slightly lower than that previously obtained by Hargreaves and Fowler (1998) (5.2 g C-CH 4 m −2 y −1 ) at Loch More. The site at Loch More bordered an extensive pool system where the water table was close to the surface, so methane emissions may genuinely have been higher at that site. However, measurements were only run there for 2 weeks and used a micrometeorological methodology, so the difference may largely be due to sampling error. The magnitude of the methane fluxes is similar to those measured at nearby sites by Macdonald et al (1996), and at the high-end of those observed in the UK (Levy et al 2012). Methane emissions at the Glencar site averaged 4.1 g C-CH 4 m −2 y −1 (Koehler et al 2011), very close to our estimate. Fluvial losses of carbon were relatively small at our site, and were also similar to those at the Glencar site (14 g C m −2 y −1 , Koehler et al 2011).

Comparison with peat core data
Combining our measurements in terms of the mass balance of carbon, we estimate a net ecosystem carbon balance −99 g C m −2 y −1 (table 1). This implies a peat accumulation rate of 2.5 mm y −1 , assuming mean values for bulk density and carbon content (0.129 g C m −3 and 52.4%, respectively from Chapman et al 2009). This result can be compared with those obtained from analysis of peat cores, where the age of peat at a given depth is measured using radiocarbon techniques, and the mass of accumulated carbon is estimated from the bulk density and carbon content of the peat. Several studies provide such data for sites within a few kilometres of the flux tower (Blackford et al 1992, Charman 1992, 1994, Coulson et al 2005, and give values averaging 20 g C m −2 y −1 (range 7-46). These values integrate over periods of several thousand years, typically 5000 y. Whether the lower values obtained by this method are because of the different time periods involved, methodological differences, or spatial sampling error, is hard to assess.
There is some uncertainty in these calculated accumulation rates, both from uncertainties in the radiocarbon dating, and in assumptions about peat bulk density and carbon content, as specific values were not always reported, and errors would propagate proportionally. As an alternative means of providing more recent dates within peat cores, the appearance of spheroidal carbonaceous particles (SCP) within peat profiles can be used to mark the onset of anthropogenic pollution from fossil fuel combustion. An analysis based on this method indicated accumulation rates of 35-209 g C m −2 y −1 at four sites in the UK (including one in northern Scotland) over the 20th century (Billett et al 2010). Our estimate based on contemporary fluxes sits in the middle of this range. We note, however, that the rate of accumulation of recent peat will not equate with the contemporary flux where there is long-term decay of the deeper peat (Charman et al 2013), and similar disagreement has been reported elsewhere (e.g. Flanagan and Syed 2011). Long-term decay of the older catotelm, may be small but not negligible in terms of calculating true accumulation rates (Clymo 1984).

Conclusions
We present a multi-year data set on the carbon fluxes at a relatively undisturbed peatland site in the Flow Country of northern Scotland. The site consistently acted as a moderate sink for CO 2 over all the measurement years. Over a range of time scales, the environmental variables which best explained the variation in NEE were quantum flux density, NDVI, temperature, and other expressions of accumulated temperature. With data of this nature, it is difficult to truly separate individual environmental effects because of co-linearity. The net uptake of CO 2 was larger than at a comparable sites at Glencar in Ireland, and Auchencorth in central Scotland, but close to the mean for other peatland and tundra sites across Europe and North America. Non-CO 2 losses comprised 11.1% of GPP, mainly from methane emissions. Considering all components of the carbon balance, the contemporary carbon sink was larger than estimates from local peat cores, based on peat accumulation over the last several thousand years. Our estimate of the contemporary flux was in the middle of the range of estimates which used SCPs to estimate peat accumulation rates over the last century.
LI-COR, Inc. in the United States and other countries. Two anonymous reviewers improved the manuscript considerably. This work was supported by the UK Natural Environment Research Council, and received funding from the European Union Seventh Framework Programme project 'Greenhouse gas management in European land-use systems', grant agreement number 244122.