Impacts of stump harvesting on carbon dioxide, methane and nitrous oxide fluxes

Abstract

During 2010-2013, we investigated the effects of stump harvesting on greenhouse gas (GHG) fluxes of carbon dioxide (CO₂), methane (CH₄) and nitrous oxide (N₂O) with the flux-gradient technique at four experimental plots in a hemiboreal forest in Sweden. All plots were clear-cut and soil scarified and two of the plots were additionally stump harvested. The two clear-cut plots served as control plots. Due to differences in topography, we had one wetter and one drier plot of each treatment. All plots exhibited substantial emissions of GHGs and we noted significant effects of wetness on CO₂, CH₄ and N₂O fluxes within treatments and significant effects of stump harvesting on CO₂ and N₂O fluxes at the dry plots. The CO₂ emissions were lower at the dry stump harvested plot than at the dry control, but when estimated emissions from the removed stumps were added, total CO₂ emissions were higher at the stump harvested plot, indicating a small enhancement of soil respiration. In addition, we noted significant emissions of N₂O at this plot. At the wet plots, CO₂ emissions were higher at the stump harvested plot, also suggesting a treatment effect but differences in wetness and vegetation cover at these plots make this effect more uncertain. At the wet plots, we noted sustained periods (weeks to months) of net N₂O uptake. During the year with simultaneous measurements of the abovementioned GHGs, GHG budgets were 1.224×10³ and 1.442×10³ gm^-2 of CO₂-equivalents at the wet and dry stump harvested plots, respectively, and 1.070×10³ and 1.696×10³ gm^-2 of CO₂-equivalents at the wet and dry control plots, respectively. CO₂ fluxes dominated GHG budgets at all plots but N₂O contributed with 17% at the dry stump harvested plot. For the full period 2010-2013, total carbon (CO₂+CH₄) budgets were 4.301×10³ and 4.114×10³ g m^-2 of CO₂-eqvivalents at the wet and dry stump harvest plots, respectively and 4.107×10³ and 5.274×10³ gm^-2 of CO₂-equivalents at the wet and dry control plots, respectively. Our results support recent studies suggesting that stump harvesting does not result in substantial increase in CO₂ emissions but uncertainties regarding GHG fluxes (especially N₂O) remain and more long-term measurements are needed before robust conclusions can be drawn.

Introduction

The terrestrial land sink absorbs about 3 Pg C annually ([10]), which contributes to limit the rate of increase in atmospheric carbon dioxide (CO₂) concentration. Forests constitute a significant part of the land sink: Tagesson et al. ([48]) used remote sensing to estimate that global forests were responsible for about 85% of the land sink or about 2.5 Pg C, which is in good agreement with an inventory estimate of 2.4 Pg C by Pan et al. ([32]). The great potential of forests to absorb CO₂ has caused increasing interests in using forests as a climate mitigation tool but forests are also important from a nature conservation point of view. This multifunctionality of forests has given forests a key role in the EU’s Green Deal ([6]), which is EU’s strategy towards sustainable development and the way forward towards net zero emissions of greenhouse gases in 2050.

How forests should be used in the climate mitigation work is however a topic of much discussion. In areas with long-lived tree species the most efficient use of forests might be to leave the forest for carbon accumulation without any management ([20]). However, in countries where forests and the use of forest products are of large importance for economy and other social benefits this is probably not a viable option. In Sweden, national strategies for reduction of greenhouse gas emissions are largely based on using residuals from forest products for substitution of fossil energy as well as using wood for replacement of cement and steel. The use of biomass for energy production in Sweden has increased continuously since the 1970s and today energy from biomass accounts for almost one third of the total energy supply. About half of this bioenergy comes from the forestry sector ([46]). The demand for bioenergy is expected to continue to increase and much of this bioenergy is expected to come from the forestry sector. It has been estimated that harvesting of logging residues could increase by a factor of 2.5 without jeopardizing sustainability ([7]).

Logging residues also include stumps, which is so far a not much utilized resource in Fennoscandia, partly because of uncertainties regarding its potentially negative environmental impacts. While no significant effects of stump harvesting on the total soil carbon (C) pools were found in the short-term to medium-term perspective ([42], [16], [33]), some studies found decreased C content in the organic layers and increased C content in mineral soil after stump harvesting ([42], [33]). This suggests a re-distribution of carbon in the soil profile as a consequence of the disturbance. It has been speculated that indirect effects (increased soil respiration due to the disturbance and mixing of the soil during stump extraction) could enhance CO₂ emissions, and thus reduce or even negate the direct effects on ecosystem respiration (lower emissions of CO₂ from the soil due to the removal of decomposable substrate), and hence reduce the potential climate mitigation benefit from using tree stumps for bioenergy production. Recent studies on soil respiration following stump harvest and/or soil scarification found no, or only small, effects of disturbance on CO₂ emissions ([40], [27], [52], [44]), although some transient effects (lower or higher emissions) were noted during the first year in some cases. It is, however, worth noticing that no effect compared to intact soil is an indication of a small indirect effect approximately equal in magnitude to the avoided stump decomposition. When scaled to plot size according to the spatial coverage of each type of disturbance, emissions from stump harvested plots were similar in magnitude to harrowed plots ([40]) and mounded plots ([41]). In contrast, Grelle et al. ([11]), using eddy covariance (EC) measurements, found a trend of increasing CO₂ emissions when comparing to soil scarified (mounded) control plots during the first three years after disturbance. During the first year, only direct effects of stump harvesting were found (i.e., lower emissions than at control plots) when the estimated (modeled) respiration of stumps left at the mounded plots were subtracted from the data. After three years, emissions were very similar at both mounded and stump harvested plots, indicating an indirect effect at the stump harvested plot and thus, a decreased climate mitigation potential of using stumps.

The full effect of stump harvesting on climate (ignoring biophysical effects, such as changes in albedo and energy exchange following clear-cutting) should not only consider the effects of soil disturbance on fluxes of CO₂ but also include fluxes of methane (CH₄) and nitrous oxide (N₂O), both powerful greenhouse gases (GHG). In general terms, soils are net sinks of atmospheric CH₄ ([39]) and the sink strength of forest soils has been estimated to be larger than for other land cover classes ([9]), which indicates that the global methane budget could be sensitive to disturbances of forest soils. Changes in hydrology following clear-cutting and stump harvesting might affect both CH₄ consumption (by methanotrophic bacteria), which takes place in aerobic parts of soils, and CH₄ production (by archaea), which might take place in anaerobic (water-saturated) parts of soils, and at anaerobic micro-sites in soils. A higher groundwater table position, as a consequence of reduced evapotranspiration following clear-cutting, may result in decreased CH₄ uptake, or even cause a site to switch from a sink to a source of CH₄ ([58], [45], [53]).

Strömgren et al. ([43]) measured CH₄ and N₂O emissions using chambers at three mesic stump harvested sites in Sweden during two years and found a dependence on GHG fluxes on disturbance type (for CH₄ and N₂O) and nitrogen (N) availability (N₂O) but fluxes were in general low. Substantial CH₄ emissions (50-70 µmol m^-2 hour^-1) were found from soil pits and wheel ruts and N₂O emissions were high (up to approx. 75 µg m^-2 hour^-1) from mounds and intact soil in areas where N availability was high. However, when scaled to plot scale no significant treatment effects were found and CH₄ and N₂O fluxes only made minor contributions (0.0-0.2% and 0.1-8%, respectively, depending on site and year) to the total GHG budgets (including CO₂). Sundqvist et al. ([45]) measured CH₄ fluxes with chambers at undisturbed, thinned, clear-cut and stump harvested plots at the ICOS Sweden station Norunda (⇒ http:/­/­www.­icos-sweden.­se/­station_norunda.­html) in central Sweden and found a reduced uptake of CH₄ in the thinned forest, as compared to the undisturbed forest, and net emissions of CH₄ at the clear-cut and stump harvested plots. The highest net emissions of CH₄ were found at the stump harvested plot.

The abovementioned studies highlight a clear need for long-term, continuous and spatially averaging measurements of GHGs at stump harvested clear-cuts in order to fill current gaps in our understanding of how stump harvesting might affect the climate. Therefore, the main objectives of this study were: (i) to study the temporally and spatially averaged dynamics of GHG fluxes following stump harvesting; (ii) to make robust inter-plot comparisons of GHG fluxes at a clear-cut and stump harvested site; and (iii) to estimate full GHG budgets following stump harvesting and relate these to “normal case scenario” in Sweden, i.e., clear-cutting and soil scarification.

Material and methods

Site description

All measurements took place at a clear-cut area (~30 ha - Fig. 1) in a hemiboreal forest near the ICOS Sweden station Norunda (60° 05′ N, 17° 29′ E) in central Sweden. Prior to harvest, the area contained a mature forest stand consisting of 58% Scots pine (Pinus sylvestris), 40.5% Norway spruce (Picea abies) and a small fraction (1.5%) of deciduous trees, mainly birch (Betula pubescens). The ground vegetation was dominated by bilberry (Vaccinum myrtillus) and feather mosses (Hylocomium splendens and Pleurozium schreberi). The soil was a highly heterogeneous sandy-loamy glacial till, rich in stones and blocks ([23]). The soils were podzolised and classified as dystric regosols (FAO Classification - [5]), with a thin (3-10 cm) organic layer.

Fig. 1 - The location of the ICOS Sweden station Norunda (60° 05′ N, 17° 29′ E) (left) and the clear-cut study site (middle). A Digital Elevation Model (DEM) showing the height above sea level for each square meter of the study area (right). The red star in the middle panel marks the position of the ICOS main tower. The black dots in the right panel represent tower positions on plots CW (lower right), CD (lower left), SHW (upper right) and SHD (upper left). Blue triangles mark the positions of 16 groundwater pipes installed down to 2.5 m depth 25 m N, E, S and W from each tower. The DEM and the aerial photograph derive from data gathered during a high spatial resolution airborne LiDAR survey in 2011. The forest was harvested in early 2009 and subsequently stump harvested, soil scarified and re-planted in May 2010.

In early 2009, the area was clear-cut and four circular experimental plots (100 m diameter) were established during the second half of 2009 and early 2010 (Fig. 1). Two of the four plots were used as control plots (i.e., clear-cut and soil scarified as this is the normal case in Sweden, see [53] for details). In May 2010, the other two plots were stump harvested with the aim to use the stump biomass for energy production. The stumps were extracted using a Steelpa head mounted on an excavator. The stumps were shaken in order to leave as much soil as possible in the pit. The Steelpa head was also used for soil scarification, creating additional mounds for planting between the planting spots created by the stump extraction. The excavator with the Steelpa head was also used for soil scarification (mounding) on the control plots where the stumps were left in the soil. All plots were planted directly following soil scarification. In total, about 2500 seedlings ha^-1 were planted (50% Norway spruce and 50% Scots pine). Changes in hydrology and fast establishment of grasses and shrubs following harvest had a negative impact on the survival of the planted seedlings. In response, additional 2200 seedlings ha^-1 (100% Norway spruce) were planted in May 2012.

Due to small differences in topography, two of the plots were wetter than the other two. In the following, “dry” and “wet” are defined based on the distance from the soil surface to the groundwater table; the closer the groundwater table was to the soil surface, the wetter it was. We had a wet control plot (CW), a dry control plot (CD), a wet stump harvested plot (SHW) and a dry stump harvested plot (SHD). This study focuses on the stump harvested plots, while the control plots will be referred to for comparisons of fluxes and GHG budgets.

Data collection and data processing

Environmental variables

A detailed description of instrumentation and location of sensors are described in Vestin et al. ([53]). Only a short outline is presented here.

A digital elevation model (DEM, cell size 1 m² - Fig. 1) was generated from a high spatial resolution airborne LiDAR survey in 2011. The DEM provided a detailed information on small-scale variations in topography at the whole clear-cut.

Profiles of soil temperature and soil moisture sensors were installed at a location considered representative for both SHW and SHD (at SHD, bordering SHW) directly following stump harvesting and soil scarification. A short mast was installed at SHD for measurements of radiation components and air temperature. Precipitation (rain and snow) was measured at 1 m height at the center of the clear-cut site.

Four groundwater pipes were installed down to 2.5 m depth at each of the four plots (Fig. 1). The groundwater level was measured continuously at one of these locations (L1) on CD from October 2010 and manually in all 16 pipes in 2010 and 2012. The differences in elevation at the ground surface between L1 and the other 15 groundwater pipes were derived from the DEM. Combined, this information made it possible to generate time series (with a daily resolution) of groundwater table positions at all pipes by assuming that fluctuations of the groundwater table position at L1 were representative for all pipes. Distance between the groundwater table and the soil surface was estimated by interpolating the groundwater surfaces for all plots using kriging in ArcGIS® ([30]). Daily average distances between the soil surface and the interpolated groundwater surfaces were calculated for each of the plots. The estimated groundwater table positions were validated using the manual measurements.

The measurements at the clear-cut were complemented with data on global radiation, incoming PAR, air temperature, air pressure, relative humidity and snow depth from the Norunda main tower site (Fig. 1), approx. 600 m from the center of the clear-cut site. When appropriate, these data were used to fill shorter gaps arising from power failures etc. in ancillary measurements at the clear-cut. In addition, data on groundwater levels measured in the mature forest in the close vicinity of the main tower were used for comparison with groundwater levels at the clear-cut plots.

Gradient measurements

Measurements of turbulence parameters and concentrations of CO₂, H₂O, and CH₄ (May 2010 - May 2013) were made at 3 m tall towers on each of the plots as detailed below. Because of instrumental problems with the original N₂O gas analyzer, N₂O measurements were only undertaken from June 2011 to September 2012.

Wind and turbulence parameters were measured at 2.5 m height by sonic anemometers (Gill Windmaster, Gill Instruments Ltd, Hampshire, UK) at a frequency of 20 Hz and were recorded on a computer. Air was continuously sampled (flow rates > 2.0 l min^-1) at two heights at each tower and pulled through 100 m long high-density polyethylene tubes (inner diameter 8 mm) to a hut placed centrally on the clear-cut. The heated air intakes had 7 µm filters and had an orifice to create underpressure in order to prevent condensation of water vapor in the tubes. The intakes were placed at 0.6 m and 2.0 m height but the lower intakes were raised to 0.8 m in June 2012 to adjust for the increasing height of the ground vegetation. The air from the intakes was analyzed for CH₄, CO₂ and H₂O with a Los Gatos tunable diode laser spectrometer (DLT-100, Los Gatos Research, CA, USA) and N₂O and H₂O (QCL Mini Monitor, Aerodyne Research, Inc., MA, USA) at a rate of 1 Hz. Each location was sampled for 15 minutes each hour, during which each level was sampled alternately for 150 s before switching to the next level. The first 45 s of data from each level were discarded and the remaining data were used to calculate an average concentration for the two levels at each tower during the 15-minute period. Hence, average concentration gradients for each of the gas species were calculated once per hour and tower. These gradients were assumed to be representative for the whole hour. Gas concentration data were quality controlled and filtered based on logical limits. Gas concentration data were also de-spiked based on standard deviation following Schmid et al. ([36]). Gas concentrations were corrected for dilution effects by water vapor. In addition, CH₄ and N₂O concentrations were also corrected for dilution by CO₂. For more details on the gradient measurements, see Vestin et al. ([53]).

Eddy covariance data

Sensible heat fluxes, friction velocities and Obukhov lengths (all derived from the sonic anemometer measurements) needed for estimating gas fluxes were calculated with the EddyPro software v. 5.0 (LI-COR, NE, USA), following standard EC methodology. Turbulent fluctuations of the raw time series data were extracted using block-averaging over 15-minute intervals matching the gradient measurements. Raw data were de-spiked following Vickers & Mahrt ([54]) before calculations of sensible heat fluxes. Data failing statistical tests for spike count, absolute limits, discontinuities and angle of attack were discarded from further analysis. The following corrections were applied: cross-wind correction for sonic temperature (applied directly by the anemometer firmware), angle of attack correction, double rotation for tilt correction, high-pass and low-pass filtering ([3]). Ogive tests showed that 15-minute averaging intervals were sufficient to capture most of the low-frequency contributions to the fluxes under turbulent conditions.

Flux-gradient calculations

The flux-gradient (FG) method, a combined eddy covariance and gradient method, was used to determine fluxes of CO₂, CH₄ and N₂O between the ecosystem and the atmosphere. The method is based on Monin-Obukhov similarity theory and is described in detail by Denmead ([8]). Vertical turbulent fluxes F_x of trace gases can be estimated as the product of the vertical concentration gradient of gas species x times its turbulent diffusivity (eqn. 1):

(1)Fx=−Kx∂Cx∂z" role="presentation"> Fx=−Kx∂Cx∂z(1) (1)Fx=−Kx∂Cx∂z

where K_x is the turbulent diffusivity (m² s^-1) of gas x, C_x is the concentration (µmol m^-3) of gas x and z is the measurement height (m). In order to determine fluxes for a specific layer of the atmosphere, eqn. 1 has to be integrated between the air sampling heights z₁ (lower) and z₂ (upper), adjusted for the displacement height d. With that, F_x becomes (eqn. 2):

(2)Fx=−κ⋅u*(Cx,z−Cx,z)ln⁡(z−dz−d)−Ψ(ς)z+Ψ(ς)z" role="presentation"> Fx=−κ⋅u*(Cx,z−Cx,z)ln(z−dz−d)−Ψ(ς)z+Ψ(ς)z(2) (2)Fx=−κ⋅u*(Cx,z2−Cx,z1)ln⁡(z2−dz1−d)−Ψ(ς)z2+Ψ(ς)z1

where κ is the von Kármán constant (set to 0.40 following [17]), u_* is the friction velocity, ψ(ς) is the integrated form of the universal function for sensible heat ([17]) and ς is a dimensionless stability parameter derived from the EC measurements. The displacement height d and roughness length z_o were estimated from wind profile measurements during the summer of 2012. The estimated d (0.2 m) and z_o (0.07 m) were assumed to be valid for all plots throughout the measurement period (2010-2013) since these parameters were mainly influenced by soil pits, mounds, piles of residues, stumps etc. and not so much by vegetation considering the poor growth of the planted seedlings. Grasses and shrubs possibly had seasonal influence on these parameters but that was not evident during the summer of 2012. The parameters were adjusted depending on the presence of snow (z_o) and on the thickness of the snow cover (d). See Vestin et al. ([53]) for more details on the FG calculations.

Footprint modeling

For the present study, it was important to select only fluxes representative for each of the plots, i.e., the according footprints should not extend outside the plot area. Horst ([13]) suggested that in the case of flux-gradient measurements, a representative flux footprint can be obtained by considering the mean height of the two gas concentration measurements. However, problems may arise if the footprints for the concentration measurements at two heights differ fundamentally at a heavily disturbed clear-cut area as was the case at the study site. Here, at a small spatial scale ( < 100 m²), the soil and surface conditions were very heterogeneous due to compaction of the soil by the harvester and other heavy machinery, piles of harvest residues and the soil scarification. In addition, turbulence parameters, which are key input variables for any footprint model, were measured only at 2.5 m height. Hence, we adopted a conservative approach and calculated flux footprint estimates for the height of the turbulence measurements. Flux footprints for smaller heights are of shorter extent, and consequently, if the modelled footprints are representative for the selected plot, it can be assumed that footprints for a smaller height would be even more representative. Running the two-dimensional footprint parameterization FFP of Kljun et al. ([18]) with input data from the four towers, hourly flux footprints were calculated for each of the plots.

Data from the airborne LiDAR survey in 2011 and from a field survey of the plot treatments were combined into a land cover classification map of the clear-cut site (Fig. 2a). From the footprint climatologies, derived by accumulating the hourly footprints for the whole measurement period, it can already be seen that the fluxes measured at each plot can be contributed to the plot’s corresponding land cover class to the largest part (up to 80%, class 1-4: plots CW, CD, SHW and SHD, respectively; class 5: clear-cut and soil scarified area surrounding all plots; class 6: harvested and soil scarified area where widely spaced seed trees were left until 2012 - Fig. 2b). Merging hourly footprints with the land cover map allowed for additional evaluation that hourly source/sink contributions per unit area to measured fluxes originated mainly (at least 80%) from the plot’s representative land cover class.

Fig. 2 - (a) A graphic representation of the land cover classifications: grey, red, blue and green fields represent footprint classes 1-4 (plots CW, CD, SHW and SHD), respectively. The yellow field represents class 5 while the dark green field represents class 6. (b) Footprint climatology, i.e., hourly flux footprints cumulated for the whole year 2012 as an example. Each of the black contour lines represents 10% contribution to the aggregated flux footprint, up to a maximum of 90% (the outermost black contour line). The red dots represent the tower locations. Due to significant overlap of the 90% footprints, only contributions up to 80% were considered. The background map depicts the tree cover derived from a LiDAR survey in 2011.

Results

Environmental conditions

Since the flux measurements started 20 May 2010, we defined a measurement year as from 20 May to 19 May in the following year. Mean air temperature were 5.0 °C, 6.9 °C and 4.9 °C for the three measurement years 2010-2011, 2011-2012, and 2012-2013, respectively. Long-term (1981-2010) average air temperature in Uppsala, approx. 25 km south of Norunda, was 6.5 °C (Uppsala University - ⇒ http:/­/­celsius.­met.­uu.­se/­). The year 2011-2012 was the warmest year during the study period, mainly due to higher air temperatures during winter months (Fig. 3a).

Fig. 3 - Environmental conditions measured at the clear-cut in Norunda from 1 June 2010 to 20 May, 2013. (a) Daily mean air temperature (°C) measured at 2 m height. (b) Daily mean soil temperature (°C) at 5 and 20 cm depth. (c) Daily sum of precipitation (mm). (d) Daily mean soil moisture (%) at 5 and 20 cm depth. e) Daily mean groundwater level at the clear-cut and at the nearby mature forest. Vertical dashed lines mark the start of a new year.

Maximum daily average soil temperatures were 19.7 °C (5 cm depth) and 17.3 °C (20 cm depth) in July 2010 (Fig. 3b), which was followed by a decreasing trend in summertime soil temperatures during the following years. Due to thick ( > 0.2 m) snow cover during the winters of 2010-2011 and 2012-2013, there was no soil frost at 5 cm depth during these years. There was no permanent snow cover on any of the plots during the winter of 2011-2012, which resulted in shallow soil frost despite the on average warmer air temperatures during this period.

Annual precipitation was 487 mm in 2010-2011, which was lower than for the following years (610 mm and 611 mm, respectively). Long-term (1981-2010) average precipitation in Uppsala was 576 mm.

Average soil moisture measured at 5 cm depth in organic soil was lower (37.8%) and more variable than in the mineral soil at 20 cm depth (45.8% - Fig. 3d). The large fluctuations in soil moisture (at 5 cm) during January-March 2011 were related to freeze-thaw events of the uppermost parts of the soil. As shown in Fig. 3e, the groundwater table at the clear-cut (measured at plot CD) was on average approximately 1 m closer to the surface than at the nearby mature forest (Norunda main station, approx. 600 m from the groundwater pipe at CD), resulting in waterlogged conditions with a large variation in space and time. As reported by Sundqvist et al. ([45]) and Vestin et al. ([53]), it is likely that the higher groundwater level at the clear-cut was a consequence of the harvest. During the latter part of 2010, it became more obvious that plots CW and SHW were substantially wetter than CD and SHD. There were small ponds scattered over all plots but large areas with superficial groundwater were found at plots CW and SHW.

The average distances from the soil surface to the interpolated groundwater surfaces were 0.16 m, 0.40 m, 0.56 m and 0.70 m at CW, CD, SHW and SHD, respectively. These findings suggest that SHW was relatively dry. However, as seen in Fig. 4, which shows the estimated average distance from the soil surface to the groundwater table for the whole site, this is likely a consequence of dry to relatively dry conditions at the slightly elevated eastern parts of SHW (cf. the DEM in Fig. 1). The groundwater table was on average close to, or above, surface at the western parts of SHW, as indicated by the blue and dark blue colors in Fig. 4. The conditions were similar at the wet control plot, CW. Note that the dark blue colors just indicate areas that potentially were flooded and not actually flooded areas. In reality, superficial groundwater often left the area as surface runoff (e.g., in the upper central part of Fig. 4, cf. Fig. 1). The results in Fig. 4highlight the highly heterogeneous hydrological conditions at the clear-cut, with scattered wet areas even at the drier plots CD and SHD.

Fig. 4 - Average distances between the soil surface and the interpolated groundwater table for the period October 2010 to May 2013 at plots CW (lower right), CD (lower left), SHW (upper right) and SHD (upper left). Dark blue colors indicate areas that could potentially be flooded but that in reality often were drained through surface runoff (which was not captured by our simplified approach with interpolated groundwater surfaces). Black circles represent the outer perimeter specified by the maximum allowed footprint extent (50 m).

Vegetation estimates from Landsat NDVI

There were no significant differences in vegetation cover between any of the plots prior to harvest in early 2009, according to the Landsat-derived summertime NDVI (Fig. 5). Directly following harvest, there were, as expected, substantial declines in NDVI at all plots, followed by a recovery phase until 2013-2014 when NDVI appeared to level off at a lower level than before the harvest at all plots. Note that summertime (JJA) NDVI data from 2013 and 2014 are taken after the flux measurement period, which ended in May 2013. The NDVI data suggest that the vegetation cover was highest at SHW, while the vegetation covers were lower but of equal magnitude at CW and SHD. The drier control plot (CD) appear to have had the lowest vegetation cover. When compared pairwise across treatments (CW vs. SHW and CD vs. SHD), the plots were significantly (p < 0.05) different in 2009 and in 2011. Pairwise comparisons within a treatment revealed significant differences (p < 0.05) in 2009-2011 (control plots) and in 2011 (stump harvested plots). The fastest initial (2009-2010) increase in NDVI were found at the wet plots CW and SHW. From 2011 and onwards, vegetation development was similar at all plots. Considering the low survival rate of the seedlings planted in 2010, most of the initial recovery in summertime NDVI was due to fast and substantial establishment of grasses and shrubs.

Fig. 5 - Landsat-derived summertime NDVI from 2005 to 2014, i.e., for both pre- and post-harvest conditions at all plots. Initially, the recovery in NDVI was faster at the wetter plots CW and SHW, indicating a faster establishment of vegetation on these plots. The vegetation development was similar at all plots from 2011 and onwards.

Net exchange of GHGs

CO₂ fluxes

When flux data were filtered based on footprint extent (i.e., non-gap-filled data where the 80% footprint extent was < 50 m) and averaged over the whole measurement period 2010 to 2013, CO₂ fluxes were 0.3 µmol m^-2s^-1 (SHW) and 0.7 µmol m^-2s^-1 (SHD). For comparison, the average CO₂ fluxes at the wet control plot CW were -0.3 µmol m^-2s^-1 and 1.0 µmol m^-2s^-1 at the dry control plot CD during the same period. The CO₂ fluxes at SHW and SHD were significantly different (p < 0.001) for the whole three-year period and for each individual year. When comparing across treatments, CO₂ fluxes at the wet plots CW and SHW were not significantly different from each other in any year, while the difference in fluxes at the dry plots CD and SHD were significant in 2011-2012 and for the period 2010-2013 as a whole.

The CO₂ fluxes were also analyzed based on wind direction (30° sectors). At SHW, average fluxes ranged between 0.0-0.6 µmol m^-2 s^-1 for footprints originating from the different wind sectors. The highest average fluxes (0.4-0.6 µmol m^-2 s^-1) originated from the most frequent wind directions (150°-240°). This sector represented a mosaic of wetter and drier conditions (Fig. 4). At SHD, average CO₂ fluxes (for each 30° wind sector) ranged between 0.4-1.5 µmol m^-2 s^-1, and were 0.4-0.7 µmol m^-2 s^-1 from footprints originating from the most frequent wind directions 150°-210°, i.e., from relatively dry areas with scattered wet spots at SHD (Fig. 4). See Vestin et al. ([53]) for details on fluxes at the control plots.

When flux data were filtered on footprint classification (i.e., gap-filled data with footprints originating from the areas specified by the land cover classification - Fig. 2a), average CO₂ fluxes at SHW and SHD were 0.7 µmol m^-2s^-1 and 0.9 µmol m^-2s^-1, respectively. As seen in Fig. 6a and Fig. 6b, there were only net emissions (daily averages) of CO₂ at SHW and SHD until the summer 2012. In contrast to the wet control plot CW, there was no daily average net uptake of CO₂ at the wet stump harvested plot SHW until the summer of 2012 (Fig. 6a and Fig. 6b). By 2012, there was daily average net uptake of CO₂ at all plots except the drier control plot CD (Fig. 6b). The cumulative NEE for the period May 2010 to May 2013 were 4.100×10³ g m^-2 of CO₂ and 4.022×10³ g m^-2 of CO₂ at SHW and SHD, respectively (Tab. 2). In comparison, the cumulative NEE at the control plots were 3.773×10³ g m^-2 of CO₂ (CW) and 5.181×10³ g m^-2 of CO₂ (CD), respectively.

Fig. 6 - Daily average GHG fluxes during 2010-2013. Black dots represent GHG fluxes at the control plots. Blue dots represent GHG fluxes at the wet stump harvested plot (SHW). Red dots represent GHG fluxes at the dry stump harvested plot SHD. Panel a (wet) and panel b (dry): CO2 fluxes ([mu]mol m-2s-1) from May 2010 to May 2013. Panel c (wet) and panel d (dry): CH4 fluxes ([mu]mol m-2hour-1) from May 2010 to May 2013. Panel e (wet) and panel f (dry): N2O fluxes ([mu]g m-2hour-1) from June 2011 to September 2012. Vertical dashed lines mark the start of a new year

Discussion

The results show large annual emissions of all studied GHGs following harvest. In contrast to the control plots, the differences in wetness did not have a large effect on total GHG emissions at the stump harvested plots but it did have impact on GPP and R_eco, as well as on the relative contribution from the different GHGs. Due to sustained periods with net N₂O uptake at the wet plot (SHW), the contribution of N₂O to the total GHG budget was reduced and was comparable to that of CH₄. In contrast, the contribution of N₂O at the dry plot (SHD) was approx. 12 times that of CH₄.

Vegetation development and its impact on greenhouse gas fluxes

Considering that fluxes of CO₂ dominated the GHG budgets at all plots (Tab. 2), the recovery time of GPP has a crucial impact on the total emissions of GHGs, both directly through photosynthesis and indirectly through increased transpiration of plants and associated changes in hydrology over time as well as through plant modulation of soil O₂ and N availability ([38]). Both cumulative GPP and NDVI data suggest more vegetation at the wet plots (CW and SHW) than at the dry plots (CD and SHD) but the relative order of the plots differ. The most negative GPP, i.e., the highest CO₂ uptake (Tab. 2), was found at plot CW for all individual years and for the period as a whole, while NDVI data (Fig. 5) indicate more vegetation at SHW throughout the measurement period. However, one cannot expect complete correlations between NDVI and GPP ([57]), and especially not for sparsely vegetated surfaces ([14]) where soil brightness can influence the NDVI signal. It is likely that stump harvesting and soil scarification disturbed the soil surface more than the soil scarification alone at the control plots ([41], [43], [44], [2]), although we did not explicitly quantify this at our plots. This might have resulted in reduced soil brightness at these plots and, as a consequence, in too high NDVI values. However, the stump harvesting could also have promoted establishment of grasses and other pioneer species, resulting in realistic estimates of NDVI at these plots.

Even though GPP was consistently more negative at the wet plots CW and SHW during 2010-2013, the recovery rates were lower than at the dry plots CD and SHD. GPP increased 2.2-fold and 3.0-fold at CW and SHW, respectively, whereas GPP increased 5.4-fold (CD) and 11.8-fold (SHD) at the dry plots. We are not aware of any studies reporting recovery rates of GPP at stump harvested clear-cuts but Humphreys et al. ([15]) reported a 2.9-fold increase in GPP three years after clear-cutting of a Douglas-fir stand in Canada. During the first three years after clear-cutting of a temperate mixed spruce and hardwood forest in northeastern USA, Williams et al. ([56]) observed a 1.7-fold increase in growing season GPP.

Greenhouse gas fluxes and treatment effects

A consequence of the very strict filtering on footprint extent (as well as on turbulence parameters, statistical tests etc.) a larger fraction of nighttime and wintertime data were removed and hence, CO₂ fluxes were biased towards uptake (CW) or lower net emissions (CD, SHD, SHW) on average. Footprints often extended outside of the 50 m radius from the towers during stable stratification and during periods with a thick snow cover (due to decreased surface roughness) and were thus removed from further analyses. On the other hand, the strict filtering, and our conservative footprint estimation approach, resulted in a strong confidence in that the observed fluxes were representative for the respective plots. By applying a paired non-parametric statistical test on data measured during the same hour, statistically significant differences were observed for fluxes of: (i) CO₂ within treatment (CW vs. CD; SHW vs. SHD) and across treatments (CD vs. SHD); (ii) CH₄ within treatment (CW vs. CD; SHW vs. SHD); and (iii) N₂O within treatment (CW vs. CD; SHW vs. SHD) and across treatments (CD vs. SHD). Thus, there were clear indications that the degree of wetness was a controlling factor of fluxes within a treatment, which highlights the importance of spatial replications of flux towers to derive robust estimates of fluxes at the ecosystem level ([31], [12]). The degree of wetness also influenced the across-treatment comparisons, indicated by the significant treatment effects of CO₂ and N₂O fluxes only at the drier plots.

By applying a less strict filtering of flux data (i.e., based on the footprint classification but with the same filtering on turbulence parameters, statistical tests etc. as described above), the surface areas of each experimental plot were maximized and total data coverage increased substantially. Through standard gap-filling and flux-partitioning procedures ([34]), it was possible to derive cumulative sums of NEE, GPP and R_eco for the period 2010-2013. We were not able to find any strong relationships between CH₄ and N₂O fluxes and typical controlling variables (e.g., soil temperature, soil moisture, groundwater table position etc.) at any of the plots. This underscores the difficulties of filling gaps in flux data at highly heterogeneous ecosystems such as a soil scarified/stump harvested clear-cut without measurements of controlling variables at microsites representing all soil surface conditions, e.g., undisturbed/disturbed soil (the latter divided into soil pits, bare mineral soil, mounds, mixed organic and mineral soil etc., stratified according to area coverage of each disturbance type), dry/wet conditions, topographical differences etc. This is especially true for N₂O fluxes that are highly variable in space and time. However, there is no general consensus of how to fill gaps in CH₄ and N₂O fluxes and the simple linear interpolation method chosen was considered sufficient to allow for robust inter-plot comparisons of these fluxes at our experimental plots.

The cumulative NEE at plots SHW and SHD were about 20% lower than fluxes reported by Grelle et al. ([11]) for a clear-cut and stump harvested Norway spruce stand in central Sweden. Grelle et al. ([11]) used continuous and spatially integrating EC measurements to investigate the effect of stump harvesting on CO₂ fluxes during 3 years. Similarly to our study, additional soil scarification (mounding) was done by the stump extractor where needed. They found a fast establishment of grasses at both stump harvested and reference plots and wet (soil moisture ranged between 30% and 69%), or even waterlogged conditions with superficial groundwater. Grelle et al. ([11]) found cumulative NEE to be in the order of 5.100×10³ and 5.220×10³ g m^-2 of CO₂ at the stump harvested plot and the mounded reference plot (i.e., clear-cut and mounded as was the case with our control plots), respectively, over the first three years after treatments. The cumulative NEE at our control plots was lower (CW, by 28%) or similar (CD). However, Grelle et al. ([11]) did not have parallel EC measurements at two plots but used annual ratios of soil respiration (measured by chambers) from stump harvested plots and reference plots to scale EC fluxes from the stump harvested plot to the reference plot.

Similarly to findings by Williams et al. ([56]) at a clear-cut White and Norway spruce forest, we noted relatively constant ecosystem respiration and a gradual recovery of GPP at all our plots. This pattern was also noted in a chronosequence study of a wide range of North American forests by Amiro et al. ([1]).

Recently, a number of studies have investigated the effects of stump harvesting and soil scarification on CO₂ fluxes in Sweden using chamber approaches. Some studies found short-term (weeks) increase in CO₂ emissions following disturbance followed by no effect, or even decreased emissions, after two years ([41], [27]). When scaled according to the area disturbed by the different treatments, Strömgren et al. ([41]) found that CO₂ emissions were slightly higher (+10%) at the stump harvested sites compared to patch scarification and similar to harrowed sites. Strömgren & Mjöfors ([40]) investigated stump harvesting compared to different methods of soil scarification in a Swedish hemiboreal forest and found no significant effect on CO₂ fluxes during the first year after treatment and lower emissions of CO₂ at disturbed sites than from intact soil during the second year. Strömgren et al. ([44]) measured CO₂ fluxes the first two years after soil scarification and stump harvesting at 14 sites in Sweden, within temperate, hemiboreal and boreal conifer forest with soils and fertility conditions common in Sweden. The main findings were significantly decreased soil respiration from both stump harvested plots and soil scarified plots compared to intact soil during the first year and lower (but not significantly lower) soil respiration during year two. No significant difference between stump harvesting and different soil preparation methods were found.

In a study at two former Norway spruce stands in Estonia, Uri et al. ([52]) reported annual emissions (extrapolated from chamber measurements during two growing seasons) of approx. 1.350×10³ and 1.800 ×10³ g m^-2 of CO₂ at stump harvested plots and approx. 1.430×10³ and 1.870×10³ g m^-2 of CO₂ at control plots, respectively, which is approx. in the range reported for the first two years in the current study (Tab. 2). In another chamber study, Webster et al. ([55]) found decreased emissions of CO₂ following biomass harvesting of a 40-year old Jack Pine stand in Ontario, Canada, with a coupling between harvest intensity and emissions; the more biomass that were removed, the lower the surface efflux of CO₂.

Combined, these aforementioned studies point towards minor effects of stump harvesting on the carbon balance, especially when compared to soil scarified clear-cuts, which is the only relevant comparison in Sweden where about 88% of the annually clear-cut area (~196.000 ha year^-1 on average 2004-2013) undergoes soil scarification ([47]). Harrowing is the most common soil scarification method used in Sweden (applied at approx. 50% of the soil scarified area, with patch scarification and mounding applied at approx. 13% each of the soil scarified area).

We are not aware of any micrometeorological studies of fluxes of CH₄ and/or N₂O at stump harvested clear-cuts and only a few chamber studies are available. Sundqvist et al. ([45]) used a system with five automatic chambers for measurements of CH₄ fluxes at undisturbed, thinned, clear-cut (at our dry control plot CD) and stump harvested (at our dry stump harvested plot SHD) sites in Norunda during short campaigns in 2010. While there was net uptake of CH₄ at the undisturbed and thinned (reduced uptake compared to the undisturbed site) sites, the chamber measurements suggested that plots CD and SHD were net sources of CH₄ (13.6 and 17 µmol m^-2 hour^-1 on average, respectively). The chambers were placed on both undisturbed and disturbed soil and fluxes ranged between -3 to 32.5 µmol m^-2 hour^-1 (CD) and -2.9 to 74.0 µmol m^-2 hour^-1 (SHD), with the highest emissions from disturbed soil. These locations coincided, in most cases, with where the groundwater table was closest to the soil surface. The results of Sundqvist et al. ([45]) were in relatively good agreement with our FG results, especially at the clear-cut control plot (Fig. 6d).

Strömgren et al. ([43]) studied the effects of soil scarification and stump harvesting on fluxes of CH₄ and N₂O using chambers during two years at three mesic sites in central Sweden. Fluxes of CH₄ were in general low, regardless of treatment, but at one of the sites significant emissions were noted from wet wheel ruts and soil pits created by the stump extraction. However, when scaled to plot level according to the spatial coverage by each disturbance type, no significant differences between treatments were found. This is in good agreement with our results where degree of wetness, rather than treatment, appear to have been the main controlling factor of CH₄ fluxes, although it cannot be ruled out that soil compaction was, at least partly, the cause of the high emissions from wheel ruts observed by Strömgren et al. ([43]).

Strömgren et al. ([43]) also found contrasting effects of disturbance type on N₂O fluxes, depending on site; highest emissions were from undisturbed soil and mounds and significantly lower fluxes from wheel ruts, mixed soil and exposed mineral soil at one of the sites while one of the other sites exhibited the highest emissions from exposed mineral soil following soil scarification. However, as in the case with CH₄, when scaled to plot scale according to the area of each disturbance type, no significant effects of either treatment on N₂O fluxes were found. This is in good agreement of our findings at the wetter plots CW and SHW, but contrasting to our findings at the drier plots CD and SHD, where we noted significantly higher emissions of N₂O (filtered for footprint extent) at the stump harvested plot. In a two-year study on the effects of site preparation of afforested grassland on fluxes of CO₂, CH₄ and N₂O from a seasonally waterlogged peaty gley soil in northeast England, Mojeremane et al. ([28]) found significantly higher CH₄ emissions after mounding, with the highest emissions from soil pits created during mounding. In contrast, mounding was found to significantly decrease N₂O emissions, with occasional net uptake of N₂O at top of mounds and in soil pits. Annual emissions ranged between 0.08×10³ and 2.63×10³ g ha^-1 of nitrogen from soil pits, mounds and undisturbed ground but no information on the spatial extent of the respective treatments are available, making it difficult to judge if the emissions were comparable to our results based on measurements that integrated over a wide range of soil surface conditions at the plots.

Similarly to Mojeremane et al. ([28]), we noted net uptake of N₂O at both wetter plots but in our case the uptake periods lasted weeks to months. The reason for this uptake is not clear but a recent study by Shurpali et al. ([38]) on agricultural land suggests that plants may affect some microbial processes connected to production and consumption of N₂O by moderating soil O₂ and N availability (both through direct plant uptake and through root exudates of carbon that stimulate heterotrophic bacteria to immobilize soil inorganic N). However, while plants may have had a role in the noted N₂O uptake at SHW during June-July 2012, it is not likely that plants had any significant impact during the periods of net N₂O uptake outside of the growing seasons. Zona et al. ([59]) reported sustained periods (weeks to a month) of N₂O uptake at a bioenergy poplar plantation at afforested agricultural land. We are not aware of any study reporting sustained periods of net N₂O uptake in hemiboreal or boreal forest ecosystems.

Greenhouse gas budgets

In order to assess the effects of stump harvesting on CO₂ fluxes at the plot level, it is desirable to measure only the effect on the heterotrophic respiration. However, when using continuous and spatially averaging flux measurements, autotrophic respiration from the vegetation is also included in the measured NEE and potential differences between plots might be explained by differences in vegetation cover and vegetation dynamics. In the following, we have used two approaches for estimating potential effects of stump harvesting on CO₂ fluxes: (i) we estimated heterotrophic respiration by subtracting autotrophic respiration from total ecosystem respiration; and (ii) we compared cumulative NEE directly.

As suggested by both CO₂ flux data and Landsat-derived NDVI, there was little vegetation at any of the plots during the first measurement year and it is reasonable to assume that heterotrophic respiration was dominating during this period. However, given the more negative cumulative GPP at the wet plots (Tab. 2), it is also reasonable to assume relatively higher contributions of autotrophic respiration mainly from grasses and shrubs at these plots and this should, if possible, be accounted for. Litton et al. ([21]) reviewed available literature and data and found that stand-level autotrophic respiration used 57% of GPP on average (range 42%-71%), across a range of forest types and was independent of stand age, resource availability and competition within a stand. We assumed that this was also valid for a clear-cut site dominated by grasses and shrubs and subtracted 57% of GPP from R_eco, the estimated heterotrophic respiration at stand level were 1.900×10³ (CW), 1.985×10³ (CD), 1.667×10³ (SHW) and 1.705×10³ g m^-2 of CO₂ (SHD) during the first year. The estimated (avoided) respiration from decomposing stumps was 0.447×10³ g m^-2 of CO₂ during the first year. The heterotrophic respiration at the wet control plot was 0.233×10³ g m^-2 of CO₂, higher (1.900×10³ - 1.667×10³ g m^-2 of CO₂) than at the wet stump harvested plot, which suggests an indirect effect (i.e., increased soil respiration as a consequence of the stump harvesting) of 0.214×10³ g m^-2 of CO₂ (0.447 ×10³ - 0.233×10³ g m^-2 of CO₂) during the first year. At the dry plots, the estimated heterotrophic respiration at CD was 0.280×10³ g m^-2 of CO₂ (1.985×10³ - 1.705×10³ g m^-2 of CO₂) higher than at SHD, indicating an indirect effect of 0.167×10³ g m^-2 of CO₂ (0.447 ×10³ - 0.280×10³ g m^-2 of CO₂) For the whole period 2010-2013, the estimated (avoided) respiration from the extracted stumps was 1.280×10³ g m^-2 of CO₂, and the net result (calculated as above) indicated indirect effects of 1.127×10³ g m^-2 of CO₂ at the wet plots and 19 g m^-2 of CO₂ at the dry plots. However, given the uncertainties in the fraction of GPP used for autotrophic respiration and the applicability of these fractions to clear-cuts dominated by grasses and shrubs, these indications of indirect effects should be interpreted with caution. In addition, GPP was estimated as the difference between measured NEE and modeled R_eco and hence, related with larger uncertainty than NEE itself.

However, also the second approach of comparing cumulative NEE points towards indirect effects. The cumulative NEE was higher at the control plots than at corresponding stump harvested plots during the first year (Tab. 2), but the differences were smaller than the estimated (avoided) respiration from the decomposing stumps. This suggests an indirect effect of 0.338×10³ g m^-2 of CO₂ at the wet plots and 0.207×10³ g m^-2 of CO₂ at the dry plots. Thus, our results for the first year are contrasting to those of Grelle et al. ([11]), who only found direct effects for the first year. During the second and third years, there was increasing influence of vegetation on the CO₂ fluxes. For the whole period 2010-2013, the net emissions of CO₂ were lower at CW (~3.773×10³ g m^-2 of CO₂) than at SHW (~4.101×10³ g m^-2 of CO₂). This difference (0.328×10³ g m^-2 of CO₂), combined with the estimated (avoided) respiration from the extracted stumps (~1.280×10³ g m^-2 of CO₂), suggests an indirect effect of stump removal at the wet plots in the order of 1.608×10³ g m^-2 of CO₂. In other words, the climate mitigation potential of these stumps was more than canceled out by increased soil respiration. However, it is possible that the difference in cumulative NEE at these plots can be explained by suppressed stump decomposition at the wetter parts of CW, while stumps decomposed “normally” at the drier parts of CW but this was compensated for by more negative GPP (45% more negative at CW than at SHW). Differences in wetness between CW and SHW (Fig. 4) also complicate the interpretation of the results. The interpolated groundwater surfaces indicated on average drier conditions at SHW than at CW. Although about 22% of the stumps remained after harvest at SHW, most of these stumps were left because of their locations close to wet areas and it is likely that decomposition of these stumps were also suppressed. The difference in cumulative NEE (1.158×10³ g m^-2 of CO₂) at the dry plots (CD and SHD) during 2010-2013 suggests an indirect effect of 0.122×10³ g m^-2 of CO₂ after three years. Considering the uncertainties in measured NEE, total harvested area, estimations of stump biomass and stump decomposition, it is not possible to determine whether or not this effect is real, but it should be viewed as an indication of reduced climate mitigation potential of tree stumps. This is in good agreement with the findings of Grelle et al. ([11]) who also found indications of reduced mitigation potential of tree stumps after three years.

During the only full year with simultaneous measurements of CO₂, CH₄ and N₂O (1 July 2011 to 30 June 2012), the total GHG budgets was dominated by CO₂ fluxes (93.5% and 81.5% contribution at SHW and SHD, respectively), with very comparable cumulative NEE at the two plots. Fluxes of CH₄ contributed little to the total budgets; only by 3.3% (SHW) and by 1.4% (SHD), while N₂O fluxes made a significant addition (17.0%) at SHD. At the wet stump harvested plot (SHW), N₂O fluxes contributed by 3.2% to the total budget. Also the carbon budgets (CO₂+CH₄ expressed as CO₂-eq.) during 2010-2013 were dominated by net CO₂ fluxes where CH₄ fluxes only contributed with 4.6% (SHW) and 2.2% (SHD).

The CH₄ fluxes (expressed as CO₂-eq.) reported here were, although relatively low, considerably higher than those reported by Strömgren et al. ([43]) from a chamber measurement study during (mainly) growing seasons at three sites in Sweden (0.0-0.2% contribution of CH₄, using GWP₁₀₀=34). The same study reported contributions of N₂O to the total GHG budgets in the order of < 1% at one site and 5-8% and 4-6% at two other sites, respectively, during the first year and 1-5% and 0.1-0.3% during the second year at these sites, respectively (using GWP₁₀₀=298). These measurements started about six months after stump harvesting and thus, the relevant comparison is fluxes during their year two and our results from 1 July 2011-30 June 2012 (i.e., roughly 13 to 25 months after stump harvest). Our results from the wet stump harvested plot (SHW) were in the same range as those reported by Strömgren et al. ([43]), while N₂O contributed substantially more at the dry stump harvested plot (SHD). Mojeremane et al. ([28]) concluded that the net effect of mounding on the total GHG budget (derived from interpolated chamber measurements), expressed as CO₂-equivalents (GWP₁₀₀ = 25 and 298 for CH₄ and N₂O, respectively) were reduced emissions the first year following mounding, and no effect in year two.

Considering that stump harvesting is not recommended at sites with moist or wet soils, our study and the above-mentioned recent studies on GHG fluxes following stump harvesting, suggest only minor indirect effects when compared to soil scarified clear-cuts. However, these minor effects mean that the climate benefit of using tree stumps for bioenergy production is likely to be reduced. No dramatic effects of disturbance on CO₂ emissions were found but there are still uncertainties regarding the importance of soil wetness, through its importance for vegetation development and related effects on GPP and evapotranspiration and through its effects on the partitioning of different GHGs in the total GHG budgets. At the drier plots, our results indicate a possible indirect effect of stump harvesting on soil respiration. In addition, this negative effect on the climate was further enhanced by substantial N₂O emissions at SHD. Also, given that stump harvesting for bioenergy production has been proposed as a measure to reduce Sweden’s CO₂ emissions in the short to medium term, one can question whether or not it is reasonable to use GWP₁₀₀ factors to evaluate the contributions from CH₄ and N₂O. By using perhaps more realistic GWP₂₀ factors, i.e., the global warming potential in a 20-year perspective (84 and 264 for CH₄ and N₂O, respectively), the relative importance of CH₄ would triple.

Conclusions

To our knowledge, this is the first study reporting continuous measurements of CO₂, CH₄ and N₂O in a clear-cut and stump harvested hemiboreal forest. It has highlighted the importance of spatial replication of flux towers and that there is a clear need for more long-term measurements that integrate GHG fluxes (especially N₂O) spatially and temporally.

Our results showed that all plots were large net sources of GHGs on an annual basis. The emissions were dominated by CO₂ and there was a clear influence on the fluxes by the establishment and dynamics of grasses and shrubs during the first years following clear-cutting and stump harvesting. Our results generally support previous studies that stump harvesting does not result in substantial increases of CO₂ emissions from soils but uncertainties remain.

The cumulative CO₂ fluxes indicated a minor indirect effect at the dry stump harvested plot, and thus, a reduced climate mitigation potential of using stumps for bioenergy production. This potential was further reduced by substantial N₂O emissions at this plot. The high N₂O emissions at the drier stump harvested plot is a cause for concern and more long-term measurements at different soils are needed.

The cumulative CO₂ fluxes at the wet plots suggested a major indirect effect at the wet stump harvested plot that led to a negative mitigation potential, i.e., the CO₂ emissions after stump harvesting were larger than the avoided respiration from the extracted stumps. However, differences in wetness, vegetation cover and distribution etc. between the wet plots make these results hard to interpret and they should be treated cautiously.

We observed a significant treatment effect on CO₂ and N₂O fluxes at the drier plots, the degree of wetness had significant effects on CO₂, CH₄ and N₂O fluxes within the same treatment. At the wet plots, we observed long periods (weeks to months) of N₂O uptake, which to our knowledge is the first such observation in hemiboreal or boreal forest ecosystems.

Acknowledgements

We are thankful for field assistance by Anders Båth and Irene Lehner at the ICOS Sweden station Norunda and by Per Weslien from the Department of Earth Sciences at University of Gothenburg. We are grateful for help with data processing by Margareta Hellström from the Department of Physical Geography and Ecosystem Science at Lund University. We thank Damien Sulla-Menashe and Eli from Department of Earth and Environment, Boston University for providing the Landsat data. We thank Dr Laura Chasmer, University of Lethbridge, for her help with processing the LiDAR data. Support for this work was provided by The Swedish Energy Agency (grant 31746-1), the Swedish Research Council (grant 621-2007-4704) and by the Linnaeus Centre LUCCI (⇒ http:/­/­www.­lucci.­lu.­se/­) funded by the Swedish Research Council. Air-borne LiDAR and digital photography for the Norunda site were acquired with support from the British Natural Environment Research Council (NERC/ARSF/FSF grant EU10-01 and NERC/GEF grant 933).

References

(1) Amiro BD, Barr AG, Barr JG, Black TA, Bracho R, Brown M, Chen J, Clark KL, Davis KJ, Desai AR, Dore S, Engel V, Fuentes JD, Goldstein AH, Goulden ML, Kolb TE, Lavigne MB, Law BE, Margolis HA, Martin T, McCaughey JH, Misson L, Montes-Helu M, Noormets A, Randerson JT, Starr G, Xiao J (2010). Ecosystem carbon dioxide fluxes after disturbance in forests of North America. Journal of Geophysical Research 115: 1-13.

(2) Andersson J, Dynesius M, Hjältén J (2017). Short-term response to stump harvesting by the ground flora in boreal clearcuts. Scandinavian Journal of Forest Research 32 (3): 239-245.

(3) Aubinet M, Vesala T, Papale D (2012). Eddy covariance - A practical guide to measurement and data analysis. Springer, Dordrecht, Netherlands, pp. 438. Online | Gscholar

(4) Axelsson E (2011). Spatiotemporal variation of carbon stocks and fluxes at a clear-cut area in central Sweden. Master Thesis, Department of Earth and Ecosystem Sciences, Lund University, Lund, Sweden, pp. 39. Online | Gscholar

(5) Chesworth W, Camps Arbestain M, Macías F, Spaargaren O, Spaargaren O, Mualem Y, MorelSeytoux HJ, Horwath WR, Almendros G, Chesworth W, Grossl PR, Sparks DL, Spaargaren O, Fairbridge RW, Singer A (2008). Classification of Soils: FAO. In: “Encyclopedia of Soil Science” (Chesworth W ed). Encyclopedia of Earth Sciences Series. Springer, Dordrecht, Netherlands, pp. 111-113. Gscholar

(6) European Commission (2019). Communication from the Commission - The European Green Deal, Brussels, Belgium, pp. 24. Gscholar

(7) De Jong J, Akselsson C, Egnell G, Löfgren S, Olsson BA (2017). Realizing the energy potential of forest biomass in Sweden - How much is environmentally sustainable? Forest Ecology and Management 383: 3-16.

(8) Denmead OT (2008). Approaches to measuring fluxes of methane and nitrous oxide between landscapes and the atmosphere. Plant and Soil 309 (1-2): 5-24.

(9) Dutaur L, Verchot LV (2007). A global inventory of the soil CH₄ sink. Global Biogeochemical Cycles 21 (4): 1-9.

(10) Friedlingstein P, Jones MW, Sullivan M, Andrew RM, Bakker DCE, Hauck J, Le Quéré C, Peters GP, Peters W, Pongratz J, Sitch S, Canadell JG, Ciais P, Jackson RB, Alin SR, Anthoni P, Bates NR, Becker M, Bellouin N, Bopp L, Chau TTT, Chevallier F, Chini LP, Cronin M, Currie KI, Decharme B, Djeutchouang L, Dou X, Evans W, Feely RA, Feng L, Gasser T, Gilfillan D, Gkritzalis T, Grassi G, Gregor L, Gruber N, Gürses O, Harris I, Houghton RA, Hurtt GC, Iida Y, Ilyina T, Luijkx IT, Jain AK, Jones SD, Kato E, Kennedy D, Klein Goldewijk K, Knauer J, Korsbakken JI, Körtzinger A, Landschützer P, Lauvset SK, Lefèvre N, Lienert S, Liu J, Marland G, McGuire PC, Melton JR, Munro DR, Nabel JEMS, Nakaoka SI, Niwa Y, Ono T, Pierrot D, Poulter B, Rehder G, Resplandy L, Robertson E, Rödenbeck C, Rosan TM, Schwinger J, Schwingshackl C, Séférian R, Sutton AJ, Sweeney C, Tanhua T, Tans PP, Tian H, Tilbrook B, Tubiello F, Van Der Werf G, Vuichard N, Wada C, Wanninkhof R, Watson A, Willis D, Wiltshire AJ, Yuan W, Yue C, Yue X, Zaehle S, Zeng J (2021). Global carbon budget 2021. Earth System Science Data Discussion 2021: 1-191. Online | Gscholar

(11) Grelle A, Strömgren M, Hyvönen R (2012). Carbon balance of a forest ecosystem after stump harvest. Scandinavian Journal of Forest Research 27 (8): 762-773.

(12) Hill T, Chocholek M, Clement R (2017). The case for increasing the statistical power of eddy covariance ecosystem studies: why, where and how? Global Change Biology 23 (6): 2154-2165.

(13) Horst TW (1999). The footprint for estimation of atmosphere-surface exchange fluxes by profile techniques. Boundary-Layer Meteorology 90 (2): 171-188.

(14) Huete AR (1988). A soil-adjusted vegetation index (SAVI). Remote Sensing of Environment 25 (3): 295-309.

(15) Humphreys ER, Andrew Black T, Morgenstern K, Li Z, Nesic Z (2005). Net ecosystem production of a Douglas-fir stand for 3 years following clearcut harvesting. Global Change Biology 11 (3): 450-464.

(16) Hyvönen R, Kaarakka L, Leppälammi-Kujansuu J, Olsson BA, Palviainen M, Vegerfors-Persson B, Helmisaari H-S (2016). Effects of stump harvesting on soil C and N stocks and vegetation 8-13 years after clear-cutting. Forest Ecology and Management 371: 23-32.

(17) Högström U (1988). Non-dimensional wind and temperature profiles in the atmospheric surface layer: a re-evaluation. Boundary-Layer Meteorology 42 (1): 55-78.

(18) Kljun N, Calanca P, Rotach MW, Schmid HP (2015). A simple two-dimensional parameterisation for Flux Footprint Prediction (FFP). Geoscientific Model Development 8 (11): 3695-3713.

(19) Lamlom SH, Savidge RA (2003). A reassessment of carbon content in wood: variation within and between 41 North American species. Biomass and Bioenergy 25 (4): 381-388.

(20) Law B, Moomaw WR (2021). Keeping trees in the ground where they are already growing is an effective low-tech way to slow climate change. The Conversation. Online | Gscholar

(21) Litton CM, Raich JW, Ryan MG (2007). Carbon allocation in forest ecosystems. Global Change Biology 13 (10): 2089-2109.

(22) Lloyd J, Taylor JA (1994). On the temperature dependence of soil respiration. Functional Ecology 8 (3): 315-323.

(23) Lundin LC, Halldin S, Lindroth A, Cienciala E, Grelle A, Hjelm P, Kellner E, Lundberg A, Molder M, Moren AS, Nord T, Seibert J, Stahli M (1999). Continuous long-term measurements of soil-plant-atmosphere variables at a forest site. Agricultural and Forest Meteorology 98 (99): 53-73.

(24) Marklund LG (1988). Biomass functions for pine, spruce and birch in Sweden. Report 45, Department of Forest Survey, Swedish University of Agricultural Sciences, Umeå, Sweden, pp. 73. Gscholar

(25) Mauder M, Foken T (2004). Documentation and instruction manual of the eddy covariance software package TK2. University of Bayreuth, Bayreuth, Bavaria, Germany, pp. 45. Gscholar

(26) Melin Y, Petersson H, Nordfjell T (2009). Decomposition of stump and root systems of Norway spruce in Sweden. A modelling approach. Forest Ecology and Management 257 (5): 1445-1451.

(27) Mjöfors K, Stromgren M, Nohrstedt HO, Gardenas AI (2015). Impact of site-preparation on soil-surface CO₂ fluxes and litter decomposition in a clear-cut in Sweden. Silva Fennica 49 (5): 1403. Online | Gscholar

(28) Mojeremane W, Rees RM, Mencuccini M (2012). The effects of site preparation practices on carbon dioxide, methane and nitrous oxide fluxes from a peaty gley soil. Forestry 85 (1): 1-15.

(29) Myhre G, Shindell D, Bréon F-M, Collins W, Fuglestvedt J, Huang J, Koch D, Lamarque J-F, Lee D, Mendoza B, Nakajima T, Robock A, Stephens G, Takemura T, Zhang H (2013). Anthropogenic and natural radiative forcing. In: “Climate Change 2013: The Physical Science Basis” (Stocker TF et al.). Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. Gscholar

(30) Oliver MA, Webster R (1990). Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Information Systems 4 (3): 313-332.

(31) Oren RAM, Hsieh C-I, Stoy P, Albertson J, McCarthy HR, Harrell P, Katul GG (2006). Estimating the uncertainty in annual net ecosystem carbon exchange: spatial variation in turbulent fluxes and sampling errors in eddy-covariance measurements. Global Change Biology 12 (5): 883-896.

(32) Pan Y, Birdsey RA, Fang J, Houghton R, Kauppi PE, Kurz WA, Phillips OL, Shvidenko A, Lewis SL, Canadell JG, Ciais P, Jackson RB, Pacala SW, McGuire AD, Piao S, Rautiainen A, Sitch S, Hayes D (2011). A large and persistent carbon sink in the world’s forests. Science 333 (6045): 988-993.

(33) Persson T, Lenoir L, Vegerfors B (2017). Long-term effects of stump harvesting and site preparation on pools and fluxes of soil carbon and nitrogen in central Sweden. Scandinavian Journal of Forest Research 32 (3): 222-229.

(34) Reichstein M, Falge E, Baldocchi D, Papale D, Aubinet M, Berbigier P, Bernhofer C, Buchmann N, Gilmanov T, Granier A, Grunwald T, Havrankova K, Ilvesniemi H, Janous D, Knohl A, Laurila T, Lohila A, Loustau D, Matteucci G, Meyers T, Miglietta F, Ourcival J-M, Pumpanen J, Rambal S, Rotenberg E, Sanz M, Tenhunen J, Seufert G, Vaccari F, Vesala T, Yakir D, Valentini R (2005). On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. Global Change Biology 11 (9): 1424-1439.

(35) Scanlon TM, Kiely G (2003). Ecosystem-scale measurements of nitrous oxide fluxes for an intensely grazed, fertilized grassland. Geophysical Research Letters 30 (16): 1-4.

(36) Schmid HP, Grimmond CSB, Cropley F, Offerle B, Su H-B (2000). Measurements of CO₂ and energy fluxes over a mixed hardwood forest in the mid-western United States. Agricultural and Forest Meteorology 103 (4): 357-374.

(37) Shorohova E, Kapitsa E, Vanha-Majamaa I (2008). Decomposition of stumps in a chronosequence after clear-felling vs. clear-felling with prescribed burning in a southern boreal forest in Finland. Forest Ecology and Management 255 (10): 3606-3612.

(38) Shurpali NJ, Rannik U, Jokinen S, Lind S, Biasi C, Mammarella I, Peltola O, Pihlatie M, Hyvönen N, Räty M, Haapanala S, Zahniser M, Virkajärvi P, Vesala T, Martikainen PJ (2016). Neglecting diurnal variations leads to uncertainties in terrestrial nitrous oxide emissions. Scientific Reports 6 (1): 410.

(39) Smith KA, Dobbie KE, Ball BC, Bakken LR, Sitaula BK, Hansen S, Brumme R, Borken W, Christensen S, Priemé A, Fowler D, Macdonald JA, Skiba U, Klemedtsson L, Kasimir-Klemedtsson A, Degórska A, Orlanski P (2000). Oxidation of atmospheric methane in Northern European soils, comparison with other ecosystems, and uncertainties in the global terrestrial sink. Global Change Biology 6 (7): 791-803.

(40) Strömgren M, Mjöfors K (2012). Soil-CO₂ flux after patch scarification, harrowing and stump harvest in a hemi-boreal forest. Scandinavian Journal of Forest Research 27 (8): 754-761.

(41) Strömgren M, Mjöfors K, Holmström B, Grelle A (2012). Soil CO₂ flux during the first years after stump harvesting in two Swedish forests. Silva Fennica 46 (1): 66. Online | Gscholar

(42) Strömgren M, Egnell G, Olsson BA (2013). Carbon stocks in four forest stands in Sweden 25 years after harvesting of slash and stumps. Forest Ecology and Management 290: 59-66.

(43) Strömgren M, Hedwall PO, Olsson BA (2016). Effects of stump harvest and site preparation on N₂O and CH₄ emissions from boreal forest soils after clear-cutting. Forest Ecology and Management 371: 15-22.

(44) Strömgren M, Mjöfors K, Olsson BA (2017). Soil-surface CO₂ flux during the first 2 years after stump harvesting and site preparation in 14 Swedish forests. Scandinavian Journal of Forest Research 32 (3): 213-221.

(45) Sundqvist E, Vestin P, Crill P, Persson T, Lindroth A (2014). Short-term effects of thinning, clear-cutting and stump harvesting on methane exchange in a boreal forest. Biogeosciences 11 (21): 6095-6105.

(46) Swedish Energy Agency (2015). Energiläget 2015 [Energy situation 2015]. Report ET 2015:08, Eskilstuna, Sweden, pp. 94. [in Swedish] Gscholar

(47) Swedish Forest Agency (2017). Statistical database on forestry. Official Statistics of Sweden, web site. Online | Gscholar

(48) Tagesson T, Schurgers G, Horion S, Ciais P, Tian F, Brandt M, Ahlström A, Wigneron J-P, Ardö J, Olin S, Fan L, Wu Z, Fensholt R (2020). Recent divergence in the contributions of tropical and boreal forests to the terrestrial carbon sink. Nature Ecology and Evolution 4 (2): 202-209.

(49) Thomas SC, Malczewski G (2007). Wood carbon content of tree species in Eastern China: interspecific variability and the importance of the volatile fraction. Journal of Environmental Management 85 (3): 659-662.

(50) Thomas SC, Martin AR (2012). Carbon content of tree tissues: a synthesis. Forests 3 (2): 332-352.

(51) Tucker CJ (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment 8 (2): 127-150.

(52) Uri V, Aosaar J, Varik M, Becker H, Kukumägi M, Ligi K, Pärn L, Kanal A (2015). Biomass resource and environmental effects of Norway spruce (Picea abies) stump harvesting: an Estonian case study. Forest Ecology and Management 335: 207-215.

(53) Vestin P, Mölder M, Kljun N, Cai Z, Hasan A, Holst J, Klemedtsson L, Lindroth A (2020). Impacts of clear-cutting of a boreal forest on carbon dioxide, methane and nitrous oxide fluxes. Forests 11 (9): 961.

(54) Vickers D, Mahrt L (1997). Quality control and flux sampling problems for tower and aircraft data. Journal of Atmospheric and Oceanic Technology 14 (3): 512-526. 2.0.CO;2" target="_blank" title="Search this paper by DOI">

(55) Webster KL, Wilson SA, Hazlett PW, Fleming RL, Morris DM (2016). Soil CO₂ efflux and net ecosystem exchange following biomass harvesting: impacts of harvest intensity, residue retention and vegetation control. Forest Ecology and Management 360: 181-194.

(56) Williams CA, Vanderhoof MK, Khomik M, Ghimire B (2014). Post-clearcut dynamics of carbon, water and energy exchanges in a midlatitude temperate, deciduous broadleaf forest environment. Global Change Biology 20 (3): 992-1007.

(57) Xiao X, Zhang Q, Braswell B, Urbanski S, Boles S, Wofsy S, Moore B, Ojima D (2004). Modeling gross primary production of temperate deciduous broadleaf forest using satellite images and climate data. Remote Sensing of Environment 91 (2): 256-270.

(58) Zerva A, Mencuccini M (2005). Short-term effects of clearfelling on soil CO₂, CH₄, and N₂O fluxes in a Sitka spruce plantation. Soil Biology and Biochemistry 37 (11): 2025-2036.

(59) Zona D, Janssens IA, Gioli B, Jungkunst HF, Serrano MC, Ceulemans R (2013). N₂O fluxes of a bio-energy poplar plantation during a two years rotation period. GCB Bioenergy 5 (5): 536-547.