Highlights

  • The bog CO2 exchange had a different spatial pattern than the underlying processes

  • Seasonal differences in CO2 exchange were linked to vegetation composition

  • Temporally changing contributions of microforms could provide stability for CO2 uptake

Introduction

Photosynthesis in boreal peatlands binds more CO2 from the atmosphere than is released back through decomposition and plant respiration. For this reason, these ecosystems are sinks of atmospheric carbon (C) and especially important as a stable storage of 273–547 Pg of C (Turunen and others 2002; Yu 2011; Kleinen and others 2012), which equals to about one-third of soil C globally. However, in certain years the annual CO2 balance may turn from uptake to release, because the C sink of boreal peatlands is small and sensitive to variations in temperature and moisture (Alm and others 1999; Waddington and Roulet 2000; Lund and others 2012) as well as light conditions (Nijp and others 2015). Due to the predicted increase in cloud cover and fluctuations in temperature and precipitation at northern latitudes (IPCC 2013), peatland C balance is likely to change. The fate of peatland C sink and storage under changing conditions has been simulated with process-based models (for example, Frolking and others 2010; St-Hilaire and others 2010; Wu and others 2011; Gong and others 2013). However, for such models more information is needed about the effect of spatially varying vegetation composition on the C cycling processes of peatland ecosystems.

Within a peatland, distance of water table (WT) to the moss surface may vary widely, resulting in different microforms—hummocks, lawns and hollows—along the WT gradient, each being habitat for a particular plant community composition. Within-site WT variation is often especially pronounced in nutrient-poor boreal peatlands, that is, bogs, dominated by dwarf-shrub vegetation at the dry end of the WT gradient, gradually changing toward sedge-dominated, wet hollows. In the moss layer, Sphagnum species adapted to certain moisture conditions dominate at different points of the WT gradient (Hayward and Clymo 1982; Rydin 1993). These microforms differ greatly in net ecosystem CO2 exchange (NEE), and one plant community may act as a C sink, whereas another is simultaneously acting as a C source (Waddington and Roulet 2000).

Both ecosystem respiration (R) by plants and decomposers and gross photosynthesis (PG) by plants have been found to be enhanced by the thicker aerobic peat layer in dry bog microforms, supporting more photosynthesizing vascular leaf area and aerobic decomposition. However, the spatial pattern of NEE is not straightforward; drier communities have been reported to have larger (Waddington and Roulet 2000), similar (Alm and others 1999; Bubier and others 2003a) and smaller (Strack and Waddington 2007) NEE as wetter ones. In previous studies, vegetation within the WT gradient is typically divided into two or three microform classes. Studying the CO2 exchange using a finer classification could reveal important thresholds for R, PG and NEE and thus help explain the contrasting spatial patterns previously observed in bogs.

Peatland plant communities differ in their phenology, and consequently, in their seasonal development of R and PG. In spring, when the deciduous sedges have low leaf area, hummocks dominated by Sphagnum mosses and dwarf shrubs already are photosynthesizing (Moore and others 2006; Riutta and others 2007). Sedges in turn have their highest photosynthesis during their midsummer peak in leaf area (Korrensalo and others 2017). The seasonal changes in R are controlled by WT (Fenner and Freeman 2011) and temperature (Lafleur and others 2005), which differ in their importance in dry and wet microforms (Maanavilja and others 2011). In addition, when moisture and temperature regimes change between years, plant communities differ in their responses; for example, in a dry year, R, PG and NEE of some plant communities may decrease, while increasing or staying similar in other communities (Bubier and others 2003b). The seasonal timing and magnitude of changes in temperature and precipitation have been found to be an important control for the interannual variation of bog C uptake (Waddington and Roulet 2000; Lund and others 2012; Peichl and others 2018). Studying the C sink of boreal bog vegetation communities over several years could reveal which seasonal patterns occur persistently and which occur as a result of interannual changes in environmental conditions, for example, only during warm or cold years. However, to our knowledge, there are only two previous studies investigating the spatial variation in bog CO2 exchange over more than two growing seasons (Strack and Waddington 2007; Pelletier and others 2011). Neither of these studies covers the whole growing season despite the potential importance of spring and autumn photosynthesis for bog CO2 uptake (Moore and others 2006; Korrensalo and others 2017).

In this study, we quantify the roles of plant communities with strongly differing species composition and WT depth for the CO2 sink of a boreal bog over three growing seasons. To fulfill our aim, we (1) quantify the response of community and ecosystem-scale NEE to variations in light, temperature, moisture and leaf area, (2) determine how plant communities differ in their seasonal timing and magnitudes of PG, R and NEE and (3) examine variations in the ecosystem-level CO2 sink by upscaling community-level NEE to an ecosystem-level estimate to be compared with ecosystem-scale CO2 flux measured by the eddy covariance method.

Methods

Study Site

The study site is an ombrotrophic bog (61°50′N, 24°12′E), located within the Siikaneva peatland complex in Ruovesi, Western Finland, in the Southern Boreal vegetation zone (Ahti and others 1968). The annual temperature sum in the area is 1318 degree days, annual rainfall is 707 mm, one-third falling as snow, and the average annual, January and July temperatures are 4.2, − 7.2 and 17.1°C, respectively (30 year averages from the Juupajoki-Hyytiälä weather station, 10 km from the site).

The large spatial variation in the WT level at the site manifests as a mosaic of vegetation communities differing in their species composition (Table 1). Dwarf shrubs (mainly Calluna vulgaris) and few small Scots pines (Pinus sylvestris) dominate the vascular vegetation in the dry end of the WT gradient. Toward the deeper WT, the dominance of field layer vegetation is shifted from dwarf shrubs to ombrotrophic sedges (Rhynchospora alba, Carex limosa) and the herb Scheuchzeria palustris. Ground-layer vegetation varies from Sphagnum fuscum dominating the high hummocks to lawns covered mainly by S. rubellum and S. papillosum and to S. majus-dominated hollows. The site also has bare peat surfaces without mosses but a sparse cover of R. alba and open water pools without any vegetation. An eddy covariance (EC) tower was located in the center of the site.

Table 1 Microform Classes Listed From the Driest to Wettest with Their Mean Sphagnum Moss, Sedge and Dwarf-Shrub Cover at the Permanent Measurement Plots and Areal Cover of the Microforms at the Study Site Based on Vegetation Inventory Within 30 m Radius from the Center
Full size table

Spatial Sampling

The variation in vegetation composition was disaggregated into seven microform classes (Table 1), open water and six vegetated microforms, mainly based on the areal cover of Sphagnum moss species, which are known to vary in relation to WT depth (Hayward and Clymo 1982; Rydin 1993). Eighteen permanent sample plots were established around the EC tower in three groups, one sample plot representing each microform (except open water) in each of the three groups. To examine how well the division into classes represents the variation in vegetation within the site, a systematic vegetation inventory was made within 30 m radius from the EC tower in July 2013. The inventory grid consisted of a total of 121 points (Appendix 1), where the percent cover of each species was visually estimated, and based on those covers, the microform class of each point was defined. The areal cover of microforms based on vegetation inventory is reported in Table 1. We used detrended correspondence analysis (DCA) to examine how plant species composition is structured and related to water table. Canonical correspondence analysis (CCA) was used to test how well the division into microform classes explains the variation in vegetation. The cover of each species was also estimated at the 18 permanent sample plots, and these data were included in both analyses to examine how well they represent the variation in vegetation within their microform class and within the whole study area. The analyses were made with Canoco 5 software (ter Braak and Šmilauer 2012).

Leaf Area Index Measurements

To follow the development of photosynthesizing leaf area at the 18 permanent sample plots over the growing seasons, we measured seasonal development of vascular green area (VGA) (Wilson and others 2007). The number of leaves of each vascular plant species in each sample plot was counted biweekly between May and September in 2012–2015. To define the VGA of each species inside the plots, the number of leaves was multiplied by the area (m2) of an average leaf on each measurement day. Average leaf size for each species was defined on each measurement day with a scanner from a sample of 10–15 leaves collected next to the plots. For Calluna vulgaris, Betula nana, Empetrum nigrum and Vaccinium oxycoccos, we measured the total length (cm) of stem containing green leaves within the plots instead of the number of leaves and average leaf area (m2) per stem centimeter instead of average leaf size.

Net Ecosystem Exchange Measurements

We used the closed chamber method (for example, Riutta and others 2007) to measure the exchange of CO2 between the vegetated land surface and the atmosphere in different environmental conditions. The measurements were conducted weekly or biweekly during the snow-free period in 2012–2014 and every three weeks in 2015. In each measurement, transparent plexiglass chamber (56 × 56 × 30 cm) was sealed airtight by placing it into the water-filled groove of an aluminum collar permanently installed in the sample plots. The CO2 concentration, photosynthetic photon flux density (PPFD) and the temperature inside the chamber were recorded every 15 s using radiation and temperature sensors and an infrared gas analyzer (EGM-4, PP systems, UK). A cooling system and a battery-operated fan maintained the chamber temperature within 2°C of the ambient temperature during the measurement. At every plot, 3–4 measurements of 90–180 s were conducted each day: first in full light, then one or two times with PPFD level reduced by 40–90% with a mosquito net, and finally with an opaque shroud. WT level was measured in perforated plastic tubes installed next to each plot, and the temperature was measured 5 and 15 cm below the moss surface. Both WT and temperature were recorded simultaneously with the other measurements at each plot.

NEE was calculated from the linear change in CO2 concentration in the chamber headspace and expressed as mg CO2 m−2 h−1. CO2 exchange measured in the dark represented the sum of autotrophic and heterotrophic respiration, R. Here, we use the sign convention where both R and PG get positive values, and positive NEE (PG minus R) indicates a sink of CO2 from the atmosphere.

Response Modeling

To reconstruct the NEE flux of the six measured microforms and to compare their responses to environmental variables, we used a nonlinear mixed-effects model with the hyperbolic light saturation curve (Larcher 2003):

$$ {\text{NEE}}_{trsij} = \frac{{Pmax_{trsi} {\text{PPFD}}_{trsij} }}{{k_{trsi} + {\text{PPFD}}_{trsij} }} - R_{trsij} + e_{trsij} $$
(1)

where NEEtrsij is the observed net CO2 exchange per square meter and the predictor PPFDtrsij is the photosynthetic photon flux density (µmol m−2 s−1) for measurement j on plot i on plot group s on Julian day r in year t. The parameters to be estimated are respiration (Rtrsi), maximum rate of light-saturated photosynthesis (\( Pmax_{trsi} \)) and the level of PPFD where half of Pmax was reached (ktrsi). The residual error (etrsi) is normally distributed with mean zero and constant variance.

In the full model, Rtrsi, \( Pmax_{trsi} \) and ktrsi were written as linear functions of fixed predictors and random effects. Possible fixed predictors were VGA, air temperature, WT depth and microform. To account for shading of vascular plants at high VGA values (Riutta and others 2007), VGA was transformed as VGAB = 1 − exp(− VGA). Air temperature was transformed as (max(T, v))2, where a second-order polynomial response was assumed for temperatures above v and no effect for temperatures below v.

The final submodels for the parameters of NEE were

$$ Pmax_{trsi} = f \, (C_{krsi} ,{\text{ VGA}}_{\text{B}} ,(\hbox{max} (T, \, v))^{2} , \, m_{t} , \, m_{kr} , \, m_{trs} , \, m_{trsi} ) $$
(2)
$$ R_{trsi} = g \, (C_{krsi} ,{\text{ VGA}}_{\text{B}} ,(\hbox{max} (T, \, v))^{2} , \, n_{t} , \, n_{tr} , \, n_{trs} , \, n_{trsi} ) $$
(3)
$$ k_{trsi} = h \, \left( {1 \, + \, a_{k} + \, a_{tr} + \, a_{trs} + \, a_{trsi} } \right) $$
(4)

where Ckrsi is the categorical predictor of microform (6 levels). T is the air temperature (°C) using a transformation v = 17 for Rtrsi and v = 12 for Pmaxtrsi. The four last terms represent the random effects. The model was fitted using the function nlme of the R program package nlme (Pinheiro and Bates 2000). A more detailed description of the response modeling is presented in Appendix 2.

Using the fixed part of the model [equations (1)–(4)], we reconstructed daily (g CO2 m−2 d−1) and seasonal (g CO2 m−2 growing season−1) cumulative PG, R and NEE fluxes at the 18 measurement plots during growing seasons 2012–2014 between May and September (Julian days 121–273), hereafter referred to as growing season. For the reconstruction, we used half-hourly PPFD and T data from SMEAR II measurement station 10 km from the site, and daily VGA of each plot calculated using the log-linear VGA development model (Wilson and others 2007). The model parameter estimates, including those of the random part, are presented in Appendix 3.

Model Validation

To validate the fixed part of our nonlinear mixed-effects model, we used the NEE measurements of growing season 2015 as an external validation dataset, which was not used for model construction. The correlation between the measured and reconstructed NEE in 2015 was explored using equations (1–4) without the random effects.

Upscaling and Comparison of Microforms

To upscale the modeled fluxes to an ecosystem-level NEE, PG and R estimates, we calculated the daily mean cumulative flux (g CO2 m−2 d−1) for each microform during the growing seasons 2012–2014 as an average of the predicted daily cumulative fluxes in the three replicate plots. Ecosystem-level daily flux was obtained as a weighted average of these microform-scale fluxes using two estimates of the contributions of microforms to the cumulative EC tower footprint as weights (Table 2). These contributions were obtained from two footprint models following Kormann and Meixner (2001) and Kljun and others (2015). Seasonal cumulative NEE, PG and R are reported in Table 2 as a mean of these two upscalings as well as 95% confidence intervals of these values based on the standard errors of the fixed part parameters (equations 1–4, Appendix 3). Open water without vegetation was found to have R close to zero (27.2 ± 23.5 mg CO2 m−2 h−1, mean ± standard deviation, n = 37) in eight small measurement campaigns conducted between May and September 2013, and thus the areal cover of water surface was included in the upscaling without a flux value by assuming that it had R and PG of zero. The small area covered by boardwalk was treated in a similar way in the upscaling. The upscalings were done for an area within about 60 m for the Kljun and others (2015) model and 130 m for the Kormann and Meixner (2001) model from the EC tower, estimated to be the radii within which 80% of the EC flux originated (see description below).

Table 2 Cumulative Ecosystem-level NEE, R/Reco and PG/GPP (g CO2 m−2 Growing Season−1) of Siikaneva Bog Site for the Growing Seasons 2012–2014 (Julian Days 121–273) Measured by Eddy Covariance Method and Upscaled from Chamber Measurements
Full size table

We compared the reconstructed seasonal R, PG and NEE (18 values of each process in each year) among years and microforms with a linear mixed-effects model with year and microforms as well as their interaction as potential fixed effects and plot group as a random effect. We used a conditional F test to test whether the fixed effects or their interaction significantly improved the model. For the fixed effects included in the model, we obtained pairwise comparisons by fitting the full model multiple times with a different year or microform acting as a default level, to which the other years and microforms were compared. This analysis was performed with function lme of the R package nlme (Pinheiro and Bates 2000).

Eddy Covariance Measurements

The upscaled daily and cumulative growing season NEE fluxes measured by the closed chamber method were compared with the NEE estimates obtained by the EC technique. The EC method offers an independent and direct estimate of the ecosystem-level exchange of energy and matter by measurement of vertical turbulent fluxes (for example, Baldocchi 2003; Aubinet and others 2012). The EC system consisted of a 3D ultrasonic anemometer (USA-1, METEK Meteorologische Messtechnik GmbH, Germany) and an enclosed H2O/CO2 gas analyzer (LI-7200, LI-COR Biosciences, USA) mounted 2.5 m above the peat surface. CO2 30-min average fluxes were obtained by post-processing the EC raw data with the EddyUH software (Mammarella and others 2016). The processing of raw data (including despiking, coordinate rotation, de-trending, time lag adjustment, spectral correction, density correction and spectroscopic correction for LI-7700) was performed following the widely accepted routines (Aubinet and others 2012). The fluxes were quality-controlled and filtered by sensor signal strength and low nighttime turbulence conditions using a threshold value of 0.1 m/s for the friction velocity (u*). According to footprint calculations with Kormann and Meixner (2001) and Kljun and others (2015) models, over 80% of the EC source area is contained within a 130 m (Kormann and Meixner 2001) and 60 m (Kljun and others 2015) radius around the EC tower in most conditions.

The NEE time series for the growing seasons 2012–2014 was gap-filled with a model (NEE = GPP − Reco). First, a respiration model was derived from the nocturnal measurements in the form

$$ R_{\text{eco}} = R_{\text{ref}} Q_{10}^{{\left( {\frac{{T_{\text{p}} - T_{\text{ref}} }}{10}} \right)}} $$
(5)

Next, the GPP was calculated as NEE-Re and modeled with a function

$$ {\text{GPP}} = \frac{{ Pmax {\text{PPFD}}}}{{k + {\text{PPFD}}}}\left( {\left( {1 - \exp \left( { - {\text{VGA}}} \right)} \right) + b} \right) $$
(6)

where the drivers include PPFD (taken from the SMEAR II station) and the annual footprint-scale VGA, with b = const describing the moss leaf area. As in the response models of chamber data, a transformation (1 − exp(− VGA)) was used to model self-shading. The resulting GPP is expressed in mg (CO2) m−2 h−1.

Results

Site Conditions

The most important gradient in species composition across the site (DCA axis 1, eigenvalue = 0.647) was related to moisture, along which the order of the plots was relative to their WT (Figure 1A, B). Vegetation formed a continuum along this gradient from high hummocks (HHU) to wet hollow (HO) and bare peat (BP) surfaces (Figure 1B). At the wet end of the moisture gradient, there was more variation also along a second gradient (Figure 1B, DCA axis 2, eigenvalue = 0.209) that separated HO and lawn (L) surfaces with almost 100% Sphagnum cover from BP without mosses. The division of plots into microform classes significantly explained the variation in species composition (CCA; pseudo-F = 32.1, p = 0.002).

Figure 1
figure 1

Mean daily air temperature (upper row), mean daytime (8AM–6PM) amount of light (photosynthetic photon flux density, PPFD, mid-row) and vascular green area (VGA) development at the six microforms (lowest row) over the three growing seasons 2012–2014.

Full size image

In all years, WT separated the microforms into three groups, and maximum growing season VGA was closely related to that division. Those groups were (1) HHU with the lowest WT and highest VGA, (2) hummocks (HU) and high lawns (HL) with intermediate WT and VGA, and (3) L, HO and BP with the highest WT (Figures 1A, 2). In the plant communities of the third group, dominating deciduous sedges made the seasonal VGA development sharper than at other communities, especially in hollows (Figure 2).

Figure 2
figure 2

(A) Mean (± standard deviation) water table (WT) relative to the moss surface during the growing seasons 2012–2014 at the six studied microforms. Negative values indicate WT below the surface. (B) DCA diagram describing the two main gradients (DCA axis 1 and 2) in species composition within the site. Samples belonging to the same microform class are marked with the same color and symbol, open symbols are used for vegetation inventory points, and solid symbols for permanent sample plots, while crosses denote species. Abbreviations of the microform classes are listed in Table 1, and species’ full names in Appendix 5.

Full size image

Year 2012 was the coldest, cloudiest and wettest with a mean growing season air T of 13.0°C, daytime PPFD of 627 µmol m−2 s−1 (Figure 2) and WT 3 cm below the moss surface (weighted average of the microform WTs in Figure 1A). Growing season 2013 was the warmest, with a mean daily temperature of 14.5°C. Water table depth and the amount of light were similar in growing seasons 2013 and 2014 (Figures 1A, 2), with WTs and daytime mean PPFDs of − 6.7 cm and 654 µmol m−2 s−1 in 2013 and − 6.5 cm and 655 µmol m−2 s−1 in 2014. In 2012–2013, daytime PPFD was the highest in the beginning and middle of the growing season, whereas in 2014, daytime PPFD showed more variation than in the previous years (Figure 2).

VGA was the highest during the warmest year (2013; Figure 2), with a site-level maximum VGA of 0.50 calculated as an average of the microform VGAs weighted by their cover within the site (Table 1). Despite the higher amount of light and deeper WT in 2014 than in 2012, VGA was rather similar in those two years (Figure 2). In 2012, the mean WT was above the surface at the three wettest microforms (Figure 1A), with low VGA at lawns and bare peat surfaces, but not at hollows in that year.

CO2 Exchange in Different Microforms

In the evaluation of the chamber model [equation (1–4)], the measured fluxes in the independent dataset of 2015 agreed well with the modeled fluxes (r2 = 0.83, slope = 0.78, intercept = − 5.5, Appendix 4). When using the model for flux reconstruction of years 2012–2014, significant differences were found among plant communities in reconstructed seasonal R (F = 1211.12, p < 0.001) and PG (F = 308.61, p < 0.001). The reconstructed R for the whole growing season decreased from the driest high hummocks to the wettest bare peat surfaces, and only lawns and hollows did not significantly differ in their seasonal R (Figure 3). There was also a larger division of plant communities into two groups, formed by the three driest and three wettest surfaces with higher and lower seasonal R, respectively (Figure 3). In reconstructed seasonal PG, the differences along the WT gradient were similar to R, except that bare peat surfaces had lower photosynthesis than the two other plant communities in the wet end of the gradient (Figure 3).

Figure 3
figure 3

Modeled cumulative growing season ecosystem respiration (R), gross photosynthesis (PG) and net ecosystem exchange (NEE) of the six vegetated microforms in growing seasons 2012–2014. The whiskers represent standard deviation among the three replicates of each microform. The microforms with different sets of letters above the bars or years with different sets of letters next to the legend symbol are significantly different (p < 0.05) according to the linear mixed-effects model. Abbreviations of the microforms as in Table 1.

Full size image

The relative differences in reconstructed seasonal NEE values among microforms were smaller than in reconstructed R and similar to reconstructed PG values (Figure 3). Significant differences were found among four groups of plant communities (F = 103.16 and p < 0.001), which were from the smallest to largest seasonal NEE: BP, HHU&HU, HL&HO and L (Figure 3). In all years, L communities were the largest CO2 sinks, with mean annual NEE varying from 194 to 287 g CO2 m−2 growing season−1 (Figure 3).

Interannual variation was higher in reconstructed PG than R (Figure 3). For all three variables—R, PG and NEE—the seasonal reconstructed fluxes were the lowest in the coldest, cloudiest and wettest year 2012 and highest in the warmest, sunniest and driest year 2013 (Figure 3). The interaction between years and plant communities was not significant in R (F = 0.58, p = 0.821), PG (F = 0.50, p = 0.877) or NEE (F = 0.49, p = 0.882).

Within-Season Variation in the Microform-Scale CO2 Exchange

For most of the growing season, plant communities kept their ranking in terms of reconstructed daily R and PG (Figure 5). This was not the case with NEE, the net of R and PG, where the absolute differences in NEE were smaller and the ranking varied over the growing seasons (Figure 5). In every spring and autumn, L communities had the highest NEE (Figure 5). In midsummer, plant communities had very similar NEE, except in the coldest and wettest year of 2012 when HL and HO had the highest midsummer NEE (Figure 5). In 2012, the different timing of maximum VGA between dry and wet communities caused the three driest plant communities to have also their seasonal NEE maximum earlier than wetter communities (Figures 1, 5). BP had much lower (mostly negative) NEE than the other communities (Figure 5).

The measured NEE fluxes followed the seasonal patterns of VGA development and air temperature (Figures 1, 4), and accordingly, the final NEE response model [equations (1–4)] included the fixed predictors of VGAB and Tair, which predicted both parameters PG and R. Of these parameters, PG had a higher threshold value v above which the temperature response is polynomial, and thus the reconstructed PG was more sensitive to Tair than reconstructed R (Figures 1, 5). This was also the case with VGAB, which had a stronger effect on PG than on R in the model (Appendix 3). As a result, R varied less seasonally than PG (Figure 5). In 2012 and 2013, the especially high midsummer VGA peak in hollows (Figure 1) resulted in their PG exceeding that of lawns in the middle of the summer (Figure 5). As PPFD is the main model driver (equation 1), short-term seasonal fluctuations in daytime PPFD were also clearly seen in the reconstructed daily PG (Figures 1, 5).

Figure 4
figure 4

Seasonal variation in net ecosystem exchange (NEE) measured in different light levels (upper row, open symbols, negative values indicate loss of CO2 from the ecosystem) and dark respiration (middle row, solid symbols) as well as medians and standard deviations of measured NEEs at each microform (lower row) in growing seasons 2012–2014.

Full size image
Figure 5
figure 5

Reconstructed daily ecosystem respiration (R), gross photosynthesis (PG) and net ecosystem exchange (NEE) at the six studied microforms over the growing seasons 2012–2014. The lines represent locally weighted scatter plot smoothing (Loess, smoothing parameter = 0.1) curves. The original variation in the modeled values is shown in Appendix 6, and the standard deviation among the three replicate plots in Figure 3.

Full size image

Interannual Differences in the Ecosystem-Scale CO2 Exchange

The ecosystem-level NEE for the whole growing season was the lowest, 71 g CO2 m−2, in the wettest and coldest year 2012 and more similar between 2013 (148 g CO2 m−2) and 2014 (114 g CO2 m−2) (Table 2). Lawn communities, with the highest areal cover within the site at 21.5%, had a high contribution to the annual NEE in all years, 55–78%. Bare peat communities had similar WT to hollows (Figure 2) but were net sources of CO2 in all years (Table 2). Both ecosystem-level R and PG were the highest in the driest year 2013, but R varied much less among the years than PG (Table 2).

Comparison of the Upscaled CO2 Exchange with Eddy Covariance Measurements

When upscaled to the ecosystem level, the daily NEE fluxes were generally smaller than NEE measured by the EC tower (Figure 6). The cumulative growing season chamber-based NEE fluxes were less than half of those measured by the EC tower in all years (Table 2). Cumulative growing season R and PG were also lower than the values derived from EC measurements (Table 2). However, the ranking in cumulative growing season NEE and PG was similar among the years in both methods so that year 2012 had the lowest and year 2013 the highest estimate (Table 2). According to both methods, there was much less interannual variability in ecosystem-scale R than in PG. The magnitude of interannual variation in NEE was smaller in upscaled chamber values (71–148 g CO2 m−2 growing season−1) than in the values based on EC measurements (173–405 g CO2 m−2 growing season−1) (Table 2).

Figure 6
figure 6

Reconstructed daily cumulative net ecosystem exchange (NEE) values upscaled to ecosystem level and NEE values measured with eddy covariance method. The lines represent locally weighted scatter plot smoothing (Loess, smoothing parameter = 0.1) curves. In 2013, there was a large gap in the EC measurements between Julian days 201 and 240, and EC flux in the latter half of 2013 is thus to a large extent based on the gap-filling model (equation 5).

Full size image

Discussion

Over the growing season, all studied microforms except BP were net sinks of CO2 from the atmosphere over the studied years, and L communities with an intermediate WT had the highest cumulative seasonal NEE, that is, CO2 sink, at our site. Similarly, in a poor fen, lawns had the highest NEE (Strack and others 2006), whereas in two earlier studies on bogs (Alm and others 1999; Moore and others 2002) the spatial differences in measured NEE were very small. However, our results contradict several other studies in different types of peatlands, where hollows have been found to be smaller C sinks than drier microforms (Waddington and Roulet 2000; Laine and others 2006, 2007; Strack and others 2006; Riutta and others 2007). On the other hand, the spatial variation in R and PG at our site gains support from previous studies; hummocks had higher R (Alm and others 1999; Laine and others 2006, 2007; Strack and others 2006) due to a thicker aerobic peat layer, and higher PG (Laine and others 2007; Munir and others 2014) due to a higher VGA than wet microforms. Many earlier studied sites have been more homogeneous, without a thick dwarf-shrub cover at the dry end of the WT gradient (Riutta and others 2007), or bogs where the hollows dominated by sedges and hollow Sphagna are either lacking (Lafleur and others 2005; Bubier and others 2003a) or replaced by open water pools (Waddington and Roulet 2000). Our results highlight the importance of heterogeneous plant species composition in regulating the spatial variability of the CO2 sink. With the fine-scale, vegetation-based microform classification used in this study, we observed two pairs of microforms, HU&HL and L&HO, that had very similar WT but distinct vegetation and different cumulative CO2 sinks in all years. In addition, in our response model microform class was a better predictor for NEE than temporally varying WT.

The high NEE of lawn communities resulted from their highest CO2 sink of all microforms in spring and autumn, a seasonal pattern in NEE that persisted over the three measured years. Lawns do not have a high VGA in spring and autumn, but they do have a very high cover of Sphagnum species adapted to intermediate WT, known to be more efficient photosynthesizers than hummocks species (Gunnarsson 2005; Granath and others 2009; Laine and others 2011). Although Sphagna are less efficient photosynthesizers than vascular plants (Moore and others 2002; Leppälä and others 2008; Otieno and others 2009; Korrensalo and others 2016), in ecosystems where VGA is low (< 1.4 m2 m−2), as in our site, mosses may account for up to 96% of the CO2 exchange (Douma and others 2007). Also, Sphagna have been found to have a strong role in the ecosystem-scale C sink in early spring, when vascular photosynthesis has not yet started (Moore and others 2006; Korrensalo and others 2017). Hollows, which have even more efficiently photosynthesizing Sphagnum species than lawns (Gunnarsson 2005; Granath and others 2009; Laine and others 2011; Korrensalo and others 2016), did not have the highest NEE in early and late summer. Hollows and lawns were rather similar in their WT, R and high Sphagnum cover. However, capitula of lawn Sphagna grow more densely than those of hollow species (Korrensalo and others 2018), which may cause the slightly higher spring and autumn PG of lawns. As observed before (Järveoja and others 2018), midsummer peak in NEE was strongly related to timing and magnitude of VGA development. Midsummer differences in NEE among plant communities were small, except for bare peat surfaces, which were mostly carbon sources.

Despite the low number of plant species, the functional diversity of bog ecosystems is high, meaning the presence of plant functional types (PFTs) with different physiology, morphology, seasonality and resource requirements (Tilman and others 1997; Cadotte and others 2008). Based on manipulation experiments (Ward and others 2009; Kuiper and others 2014) or observed natural variation (Alm and others 1999; Bubier and others 2003b), bog plant communities and PFTs are known to have differential responses to changes in environmental conditions. The functional diversity of bogs has been found to make their C sink more stable over a growing season than the C sink of fens, which often have functionally more homogeneous, sedge-dominated vegetation (Bubier and others 1998; Leppälä and others 2008). Different responses of microforms with differing vegetation were also observed in this study, where contributions of plant communities to the ecosystem-scale CO2 sink varied both within the growing seasons and between years with different temperature, WT and PPFD. Bog PFTs are known to respond differently to changing environmental conditions. Sphagna growing in the wet microforms respond with changed growth rate both to WT and to temperature variations more readily than hummock species (Robroek and others 2007). The abundance and biomass production of vascular PFTs are also known to react in opposite directions as a result of WT drawdown (Mäkiranta and others 2018), as deciduous vascular plants have generally less transpiration-decreasing morphological adaptations than evergreen dwarf shrubs. Over wide environmental gradients, bog species composition may change completely, but at the same time, the functional identity of bog ecosystems may remain rather similar (Robroek and others 2017). Nonetheless, more research is needed to determine whether the large diversity of PFTs and the varying responses of vegetation communities could, over timescales from several years to decades, make the C sink of bogs more stable than of peatlands with more functionally homogeneous vegetation composition, and which mechanisms would be critical for this stability. Studies reporting C fluxes for multiple years in varying types of peatlands (Peichl and others 2014; Helfter and others 2015; Lund and others 2015) could enable a multiyear comparison of peatland ecosystems of similar climatic conditions but different vegetation structure. Further, the variability of plant species and functional traits in bog ecosystems should be studied in relation to their C sink properties.

The site-level NEE measured with either EC or chamber methods fell within the range of several previously measured boreal bogs (Laine and others 2007; Koehler and others 2011; Petrescu and others 2015; Wilson and others 2016). It was lower than measured in a temperate bog with more shrubs and higher levels of PPFD (Lafleur and others 2001, 2003), and higher than the NEE of 67.8–72.2 g CO2 m−2 a−1 measured in a coastal blanket bog (Lund and others 2015). Because our site was rather wet for a bog, a deeper summertime WT increased the growing season CO2 sink measured by EC tower from 173 to 405 g CO2 m−2 growing season−1 between years 2012 and 2013, instead of decreasing it, as in drier bogs (Roulet and others 2007; Helfter and others 2015).

The different magnitude of fluxes estimated by chambers and EC tower raises questions about the possible error sources of the two methods. A mismatch between upscaled chamber fluxes and EC tower flux, although smaller than ours, was also found by Maanavilja and others (2011), where the high NEE peaks observed by the EC tower were not seen in the upscaled chamber fluxes. In chamber data, the choices made in modeling may cause different outcomes (Laine and others 2009). In our NEE model, WT was only included as a predictor in the form of microform, although WT is known to be an essential environmental factor explaining R, PG and NEE in peatlands (for example, Tuittila and others 2004; Laine and others 2006; Riutta and others 2007). However, the seasonal rhythm of WT variation did not differ among plant communities (p > 0.05 in all years, data not shown), and microform also includes the variation related to species composition. Additionally, our NEE model is supported by the model validation. In EC measurements, comparison of several towers at the same site has revealed that spatial variability causes substantial differences in the measured CH4 flux magnitude and diel course (Peltola and others 2015), but as far as we know such analysis has not been conducted on CO2 fluxes. Gap filling is an important source of uncertainty in cumulative EC fluxes (Moffat and others 2007), and most of this uncertainty originates from the different treatments of long gaps, which in our case occur mainly at the beginning and end of the growing season. Based on the uncertainty of the EC gap filling method used in this study, reported in Moffat and others (2007), the 95% confidence intervals for EC-based seasonal NEE values can be estimated to be approximately 137–209, 348–462 and 318–404 for growing seasons 2012, 2013 and 2014, respectively. In addition, measured EC fluxes are prone to systematic and random uncertainties, usually ranging between 10 and 30% of the cumulative seasonal flux (for example, Baldocchi 2003; Rannik and others 2006, 2016). Considering the confidence intervals of EC-based and chamber-based NEE (Table 2), the uncertainty ranges of the seasonal values based on the two methods did not overlap. The difference between EC-based NEE and chamber-based NEE was the largest in 2013, when it was higher than the mean annual CO2 sink reported for northern (arctic to temperate) open peatlands (Petrescu and others 2015). Further, in our data interannual variability in NEE was larger for EC data (173–405 g CO2 m−2 growing season−1) than for upscaled chamber fluxes (77–141 g CO2 m−2 growing season−1). The difference between the two methods was thus substantial for interpreting both the magnitude and the interannual patterns of NEE. Based on these results, the observed difference between the most widely used methods for NEE estimation warrants more research to reliably predict C sink of northern peatlands. This research should include more work on revealing the magnitude of error sources of the CO2 flux measurements and data processing by both EC and chamber methods. In this process, these two methods will provide an invaluable comparison against each other’s.

Conclusions

Microforms along the WT gradient of an ombrotrophic boreal bog had strongly varying species composition and large variation in their seasonal reconstructed PG and R. Species composition has the potential to regulate CO2 sink strength of the microforms. The results show for the first time that two microforms with similar WT but distinct vegetation had a different seasonal NEE for three consecutive years. The seasonal reconstructed NEE was the highest in the most abundant microform, lawn. The seasonally and interannually varying contributions to ecosystem-scale NEE of microforms dominated by different plant functional types suggests this variability in vegetation composition may add to the temporal stability of the site CO2 sink.

The reconstructed daily chamber fluxes were less than half of the EC fluxes, and the interannual variability in seasonal ecosystem-level NEE differed between the two methods. This result calls for exploring more the uncertainties related to both of the methods.