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.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

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):

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):

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.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Flux data filtering

Data was lost due to snow and frost on the sonic anemometers during wintertime, due to heavy rain events, power failures, calibrations and other instrument maintenance. In total, 15% (CW), 24% (CD) 12.5% (SHD) and 19.0% (SHW) of the data were lost for these reasons.

The Mauder & Foken ([25]) data quality tests on developed turbulent conditions and steady state conditions were applied, resulting in a quality flag system (0-1-2) of the sensible heat flux data. Data with flag 2 were discarded from further analyses.

Total errors of the GHG fluxes were estimated as the quadratic sums of errors in fluxes caused by the gradient estimations and by errors in fluxes caused by estimations of the turbulent diffusivities. Data were discarded when the total error of the GHG fluxes was larger than 2 µmol m^-2 s^-1 (CO₂), 20 µmol m^-2 hour^-1 (CH₄) or 40 µg m^-2 hour^-1 (N₂O).

The dependence of fluxes on friction velocity was evaluated with the R (version 3.3.2) package REddyProc (version 0.8-2) and were found to vary slightly between plots and years, due to differences in instrumentation (sonic anemometers) and physical properties of each plot. Average u_* thresholds of 0.16 ms^-1, 0.15 ms^-1, 0.14 ms^-1 and 0.13 ms^-1 where used for filtering out low-turbulence data at plots CW, CD, SHW and SHD, respectively, during 2010-2013.

The GHG flux data were further filtered according to the flux footprints in two separate ways: footprint extent and footprint classification. Firstly, hourly flux data were filtered according to the 80% footprint extent for insuring a relatively homogeneous source area, i.e., fluxes with an 80% footprint extending outside the 50 m radius were discarded from further analyses. This approach minimized any possible contributions to the measured fluxes from other land-cover classes and allowed for a more direct comparison of net GHG fluxes at each plot, at the cost of total data coverage. This filtering approach is in the following referred to as footprint extent. Secondly, flux data were filtered according to the footprint land cover classification, which, since it was based on the real plot geometry, maximized the plot areas (see Fig. 2b and Section “Footprint modeling” for details). Small contributions (<10%) to the measured fluxes that originated from other land cover classes were allowed. This approach minimized data loss and these data were then gap-filled to derive GHG budgets. This filtering approach is in the following referred to as footprint classification. Despite the very strict filtering on footprint extent (and on all other criteria mentioned above), 21.3-30.6% of data across all GHGs and both plots still remained for direct inter-plot comparisons. The data filtering on footprint classification resulted in 37.3-57.3% of data remaining for creating GHG budgets on both plots. Possible problems with temporal autocorrelation of the fluxes were reduced by selecting data (filtered for footprint extent) measured during the same hour on both plots. A non-parametric paired statistical test (Wilcoxon signed rank test for zero median) for evaluation of possible treatment effects on GHG fluxes was carried out using Matlab.

All gradients, turbulent diffusivities, statistics and final GHG fluxes were calculated with Matlab R2015b (Mathworks, Inc., MA, USA).

Gap-filling, flux-partitioning and global warming potential

CO₂ flux data

Gaps in CO₂ flux data were filled using the R software v. 3.3.2 and the R package REddyProc v. 0.8-2. REddyProc is a gap-filling and flux-partitioning tool that applies the methods of Reichstein et al. ([34]). In short, gaps in data are filled with average values during similar meteorological conditions (or during the same time of the day in case no meteorological data were available) for a time window of varying size.

Only original data passing quality tests and filtering were used for flux-partitioning. The flux-partitioning assumes zero gross primary productivity (GPP) during night-time, i.e., the measured net ecosystem exchange (NEE) equals total ecosystem respiration (R_eco) when global radiation < 20 Wm^-2. The Lloyd & Taylor ([22]) model was then used to derive an empirical model of R_eco based on air temperature. Ecosystem respiration was extrapolated to daytime conditions using daytime air temperature. GPP was estimated as the difference between measured NEE and modeled R_eco. See Reichstein et al. ([34]) for a detailed description of the gap-filling and flux-partitioning procedures.

CH₄ and N₂O flux data

Gaps in CH₄ were filled using linear interpolation of daily average methane fluxes. Daily average methane fluxes were calculated for days with at least six hourly measurements passing all filtering criteria and where outliers (σ>5) were removed prior to the averaging. Gaps in N₂O flux data were filled using linear interpolation of daily averaged data, similar to Scanlon & Kiely ([35]). Daily averages were calculated in the same way as for the CH₄ data.

Global warming potential

Gap-filled CH₄ and N₂O fluxes were also converted into CO₂-equivalents by multiplication with a global warming potential factor in a 100-year perspective (GWP₁₀₀) of 28 (CH₄) and 265 (N₂O) ([29]). These GWP₁₀₀ factors do not include any climate-carbon feedbacks (see [29] for definitions) and hence, represent a conservative climate forcing estimate.

Vegetation development

The vegetation cover on each of the experimental plots was estimated during the main growing seasons (June, July and August) of 2005 to 2014 by using Landsat images with a spatial resolution of 30 × 30 m. The growing season Normalized Difference Vegetation Index (NDVI - [51]) of each plot was calculated by averaging all high-quality Landsat data for 3 × 3 pixels during summer (JJA) of each year of the study period. These 90 × 90 m areas overlapped the areas specified by the footprint classification (Fig. 2a) by 70-85%. See Vestin et al. ([53]) for more details.

Estimations of stump biomass and stump decomposition

A part of the original forest stand, directly to the East of the experimental plots (Fig. 2a, Fig. 2b), was not harvested at the same time as the experimental plots were harvested (early 2009). This made it possible to derive linear relationships between stem diameters at breast height (DBH, ~1.3 m height) and stump diameters D at 20, 30 and 40 cm height for Norway spruce and Scots pine in this part of the stand (Tab. 1). In total, 16 trees of each species were measured (mean cross-calipered diameters) during late 2009 (see [4] for more details on these measurements). The derived linear relationships were assumed to be valid for the pre-harvest forest on plots SHW and SHD, allowing a derivation of a theoretical DBH for the stumps. All stumps on plots SHW and SHD were inventoried (species, mean height and cross-calipered mean diameter) before and after the stump harvest in May 2010, which, in combination with biomass functions for the stumps including coarse roots (Scots pine: e^{11.1106(DBH/(DBH+12))-3.3913}; Norway spruce: e^{10.5381(DBH/(DBH+14)-2.4447} - [24]), made it possible to estimate how much biomass (dry biomass and carbon only) was removed from SHW and SHD when the stumps were extracted.

Lamlom & Savidge ([19]) reported that the commonly assumed carbon content of woody biomass (50%) is an over-simplification and suggested values in the range of 47.2% to 55.2% for North American conifers, with an average of 51.1%. Thomas & Malczewski ([49]) reported 50.8% as an average value for conifers in Eastern China and a literature review by Thomas & Martin ([50]) suggested 50.8 ± 0.6% as an average value for conifers in temperate and boreal regions. However, due to a lack of species-specific values for Scots pine and Norway spruce, we assumed 50% carbon content of the dry woody biomass.

Melin et al. ([26]) studied decomposition of 0-39 years old Norway spruce stump and coarse root samples from three locations along a latitudinal transect in Sweden. Their results suggested a relative decomposition rate of 4.6% year^-1 of the woody biomass in stumps and coarse roots and that there was no need to include any time lag to account for the time it might take for microbes to colonize the stump-root system (i.e., it is reasonable to assume that stump-root decomposition starts the same year the stems were harvested). Shorohova et al. ([37]) reported relative decompositions rates of Norway spruce and Scots pine to be 5.2% year^-1 and 4.8% year^-1, respectively, in Southern Finland. However, since this study did not include fresh stumps (only 5-40-year-old stumps), and did not include root systems, we applied the relative decomposition rate suggested by Melin et al. ([26]) to both spruce and pine stumps in our study to estimate how much the extracted stumps would have contributed to the CO₂ budgets at SHW and SHD if they had been left to decompose in the soil.

  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.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

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).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

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.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

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): CO₂ fluxes (µmol m^-2s^-1) from May 2010 to May 2013. Panel c (wet) and panel d (dry): CH₄ fluxes (µmol m^-2hour^-1) from May 2010 to May 2013. Panel e (wet) and panel f (dry): N₂O fluxes (µg m^-2hour^-1) from June 2011 to September 2012. Vertical dashed lines mark the start of a new year

  Enlarge/Shrink   Download   Full Width  Open in Viewer

While the cumulative NEE was similar at SHW and SHD, there were large differences in how the fluxes were partitioned into GPP and R_eco. During 2010-2013, SHW had higher ecosystem respiration than SHD (6.594×10³ g m^-2 of CO₂ vs. 5.392×10³ g m^-2 of CO₂) as well as a more negative GPP (-2.493 ×10³ g m^-2 of CO₂ vs. -1.370×10³ g m^-2 of CO₂). At the control plots, R_eco were 7.385×10³ g m^-2 of CO₂ (CW) and 6.792×10³ g m^-2 of CO₂ (CD), while GPP were -3.612×10³ g m^-2 of CO₂ (CW) and -1.611×10³ g m^-2 of CO₂ (CD). See Tab. 2for details during individual years.

CH₄ fluxes

Average CH₄ fluxes (filtered on footprint extent) for the whole period 2010-2013 were 21.1 µmol m^-2 hour^-1 and 7.2 µmol m^-2 hour^-1 at SHW and SHD, respectively. Corresponding CH₄ fluxes at the control plots were 26.1 µmol m^-2 hour^-1 (CW) and 8.2 µmol m^-2 hour^-1 (CD). The CH₄ fluxes were significantly different (p<0.001) for the whole period 2010-2013 and for each individual year. When compared to the control plots (CW vs. SHW and CD vs. SHD), fluxes were not significantly different in any of the individual years and not for the period as a whole. When averaged over 30° wind sectors, CH₄ fluxes ranged between 12.5-34.6 µmol m^-2 hour^-1 at SHW, with the lowest net emissions (12.5-18.3 µmol m^-2 hour^-1) in the sector 90°-270°. This sector encompassed both some of the driest parts of the clear-cut and some areas with superficial groundwater (Fig. 4). Fluxes from the other 30° sectors ranged between 24.3-34.6 µmol m^-2 hour^-1. There were smaller differences between sectors at the dry stump harvested plot SHD, where average fluxes ranged between 3.7-8.6 µmol m^-2 hour^-1. The net emissions were slightly higher from the sector 120°-210°, which, although relatively dry on average, also contained some of the wettest parts at SHD (Fig. 4).

The average gap-filled CH₄ fluxes (i.e., data filtered on footprint classification) were 16.8 µmol m^-2 hour^-1 and 6.6 µmol m^-2 hour^-1 at SHW and SHD, respectively for the period 2010-2013 (Fig. 6c and Fig. 6d). Corresponding CH₄ fluxes at the control plots were 22.6 µmol m^-2 hour^-1 (CW) and 8.2 µmol m^-2 hour^-1 (CD - Fig. 6c and Fig. 6d). Directly following stump harvesting and soil scarification in May 2010, daily average CH₄ fluxes at SHW were high (up to approx. 100 µmol m^-2 hour^-1) and variable, with emission peaks that correlated with precipitations events and increased soil moisture levels in the topsoil (Fig. 3). The net emissions decreased during the fall 2010 and remained relatively low, with the exception of smaller emission peaks related to fluctuations in groundwater table position in April and in June 2011, until the summer of 2012. Initially, there was net consumption of CH₄ at SHD in the order of -20 µmol m^-2 hour^-1 but this turned into net emissions in mid-July 2010. SHD exhibited low and relatively steady net emissions of CH₄ throughout the rest of the measurement period, but there were signs of small emission peaks related to changes in hydrology (e.g., a small emission peak at the end of March 2011 during snow melt). See Vestin et al. ([53]) for details on fluxes at the control plots.

In terms of the total carbon budget for 2010-2013, the cumulative CH₄ fluxes at SHW and SHD were negligible and only contributed with 5.4 g m^-2 of carbon and 2.5 g m^-2 of carbon, respectively. In comparison, the control plots contributed with 8.9 g m^-2 of carbon (CW) and 2.5 gm^-2 of carbon (CD).

N₂O fluxes

Average N₂O fluxes (filtered for footprint extent) during the period June 2011 to September 2012 were 31.6 µg m^-2 hour^-1 (SHW) and 111.0 µg m^-2 hour^-1 (SHD), while the N₂O fluxes at the control plots were 47.4 µg m^-2 hour^-1 (CW) and 64.9 µg m^-2 hour^-1 (CD). The fluxes were significantly different (p<0.001) for the whole measurement period. When compared across treatments (CW vs. SHW and CD vs. SHD), differences were only significant (p<0.001) for the drier plots CD and SHD. Average fluxes per 30° wind sector at SHW ranged between 17.4-38.8 µg m^-2 hour^-1, with the lowest fluxes for footprints originating from the relatively dry sector 60°-90° and the relatively wet sector 270°-330° (Fig. 4). At SHD, there were no clear dependence on wind direction for the average N₂O fluxes, which ranged between 103.6-121.5 µg m^-2 hour^-1.

The N₂O fluxes filtered for footprint classification and averaged over the whole measurement period June 2011 to September 2012 were 20.6 µg m^-2 hour^-1 (SHW) and 101.0 µg m^-2 hour^-1 (SHD - Fig. 6e and Fig. 6f). Corresponding N₂O fluxes at the control plots were 37.4 µg m^-2 hour^-1 (CW) and 57.5 µg m^-2 hour^-1 (CD - Fig. 6e and Fig. 6f).

There was a decreasing trend in N₂O fluxes at both plots in June-July 2011, that coincided with a decreasing trend in soil moisture (5 cm) and groundwater level. This was turned into emissions peaks at the end of July (SHW) and in the beginning of August (SHD) following precipitations events that affected both soil moisture (5 cm) and groundwater level (Fig. 3). Daily average N₂O fluxes at SHD remained high, reaching approx. 200 µg m^-2 hour^-1 at several occasions, and with large variability throughout the measurement period. Similar to the wet control (CW), daily average N₂O fluxes at SHW were lower and less variable than at the corresponding dry plots. During November 2011, there was a switch from net N₂O emissions into substantial net N₂O consumption (approx. -40 µg m^-2 hour^-1) at SHW. Fluxes were low and switched between net consumption and net emission from the end of December 2011 when the soil repeatedly froze and thawed (at 5 cm depth). The period with net consumption, or close-to-zero fluxes, of N₂O lasted until March 2012 when it was turned into a small emission peak correlated with freeze-thawing events (Fig. 3) of the soil at 5 cm depth. This was followed by other periods with net consumption of N₂O (April, May and June-July 2012). See Vestin et al. ([53]) for details on fluxes at the control plots.

Cumulative fluxes for a full year (defined as 1 July 2011 to 30 June 2012) were 0.15 g m^-2 of N₂O (SHW) and 0.93 g m^-2 of N₂O (SHD). This translates into annual emissions of 1.0×10³ and 5.1×10³ g ha^-1 of nitrogen at SHW and SHD, respectively. Corresponding values for the control plots were 0.29 g N₂O m^-2 (CW) and 0.49 g N₂O m^-2 (CD), which equals 1.9×10³ g ha^-1 (CW) and 3.1×10³ g ha^-1 (CD) of nitrogen. The cumulative fluxes for the whole period June 2011 to September 2012 were 0.22 g m^-2 of N₂O (SHW) and 1.11 g m^-2 of N₂O (SHD), which equals 1.4×10³ and 7.0×10³ g ha^-1 of nitrogen, respectively.

Total GHG budgets

Gap-filled CH₄ and N₂O fluxes were summed over the whole measurement period and converted into grams of CO₂-equivalents per square meter in order to derive full GHG budgets for all plots and to compare the relative importance of the different GHGs at a specific plot. Due to technical problems with our initial N₂O analyzer, it was not possible to derive full GHG budgets for the full three-year period. In Tab. 2, we present a complete GHG budget for a full year, defined as from 1 July 2011 to 30 June 2012 (marked *** in Tab. 2) in order to facilitate comparisons of the relative importance of the different GHGs. The total GHG emissions, in a 100-year perspective, were 1.224×10³ g m^-2 of CO₂-eq. And 1.442 ×10³ g m^-2 of CO₂-eq. at SHW and SHD, respectively. Corresponding budgets at the control plots were 1.070×10³ g m^-2 of CO₂-eq. (CW) and 1.696×10³ g m^-2 of CO₂-eq. (CD). For the full measurement period 2010-2013, the cumulative net emissions of CO₂ and CH₄ were 4.301×10³ g m^-2 of CO₂-eq. (SHW) and 4.114×10³ g m^-2 of CO₂-eq. (SHD), while the combined CO₂ and CH₄ fluxes at the control plots were 4.107×10³ gCO₂-eq. m^-2 (CW) and 5.274×10³ gCO₂-eq. m^-2 (CD). Despite the much shorter measurement period, N₂O fluxes still made significant contributions to the GHG budget, both in individual years and for the whole period 2010-2013 (6.7%) at the SHD (see total sums marked ** in Tab. 2).

Stump biomass and stump decomposition

The stump inventory showed that the pre-harvest forest stand density was approx. 730 trees ha^-1. In total, there were 681 Scots pine, 474 Norway spruce and 16 Birch at plots SHW and SHD combined (a total area of approx. 1.6 ha was stump harvested). After harvest, 22.1% of the initial stumps remained (25.4% of Scots pine, 17.7% of Norway spruce and 12.5% of Birch).

In total, 84.9 Mg of stump biomass was harvested (53.1 Mg ha^-1), which translates into average value of 2.65×10³ g m^-2 of carbon for the stump harvested plots (Tab. 3). The estimated contribution to R_eco from the stumps, if they had been left in the soil, was 1.280×10³ g m^-2 of CO₂ for the three-year period (Tab. 3).

  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