Impacts of Clear-Cutting of a Boreal Forest on Carbon Dioxide, Methane and Nitrous Oxide Fluxes

The 2015 Paris Agreement encourages stakeholders to implement sustainable forest management policies to mitigate anthropogenic emissions of greenhouse gases (GHG). The net effects of forest management on the climate and the environment are, however, still not completely understood, partially as a result of a lack of long-term measurements of GHG fluxes in managed forests. During the period 2010–2013, we simultaneously measured carbon dioxide (CO2), methane (CH4) and nitrous oxide (N2O) fluxes using the flux-gradient technique at two clear-cut plots of different degrees of wetness, located in central Sweden. The measurements started approx. one year after clear-cutting, directly following soil scarification and planting. The study focused on robust inter-plot comparisons, spatial and temporal dynamics of GHG fluxes, and the determination of the global warming potential of a clear-cut boreal forest. The clear-cutting resulted in significant emissions of GHGs at both the wet and the dry plot. The degree of wetness determined, directly or indirectly, the relative contribution of each GHG to the total budgets. Faster establishment of vegetation on the wet plot reduced total emissions of CO2 as compared to the dry plot but this was partially offset by higher CH4 emissions. Waterlogging following clear-cutting likely caused both plots to switch from sinks to sources of CH4. In addition, there were periods with N2O uptake at the wet plot, although both plots were net sources of N2O on an annual basis. We observed clear diel patters in CO2, CH4 and N2O fluxes during the growing season at both plots, with the exception of CH4 at the dry plot. The total three-year carbon budgets were 4107 gCO2-equivalent m−2 and 5274 gCO2-equivalent m−2 at the wet and the dry plots, respectively. CO2 contributed 91.8% to the total carbon budget at the wet plot and 98.2% at the dry plot. For the only full year with N2O measurements, the total GHG budgets were 1069.9 gCO2-eqvivalents m−2 and 1695.7 gCO2-eqvivalents m−2 at the wet and dry plot, respectively. At the wet plot, CH4 contributed 3.7%, while N2O contributed 7.3%. At the dry plot, CH4 and N2O contributed 1.5% and 7.6%, respectively. Our results emphasize the importance of considering the effects of the three GHGs on the climate for any forest management policy aiming at enhancing the mitigation potential of forests.


Introduction
Forests play an important role in the mitigation of fossil fuel emissions by sequestering atmospheric carbon dioxide. Pan, Birdsey [1] estimated the global annual forest sink to be 2.4 ± 0.4 Pg C on average 575. 6 mm (1981-2010). The soil is a sandy-loamy glacial till with a high content of stones and blocks [72], covered with a thin organic layer (3-10 cm thick). The pre-harvest forest was a mixed Scots pine (Pinus sylvestris, 58%) and Norway spruce (Picea abies, 40.5%) stand, with a small fraction of birch (Betula sp., 1.5%). Before harvest, the ground vegetation was dominated by bilberry (Vaccinum myrtillus) and feather mosses (Hylocomium splendens and Pleurozium schreberi). The forest stand (~30 ha) was harvested in early 2009 ( Figure 1) and was subsequently soil scarified (mounding or harrowing) and replanted in May 2010. Approximately 2500 seedlings ha −1 were planted (50% each of Scots pine and Norway spruce). Following harvest, grasses and shrubs quickly dominated the site, which, together with changes in hydrology, had a negative impact on the survival rate of the seedlings. Hence, in May 2012, an additional 2200 seedlings ha −1 (100% Norway spruce) were planted.

Environmental Variables
Directly following soil scarification in May 2010, sensors for measurements of soil temperature and soil moisture were installed at different depths at a location considered to be representative for both plots (at plot 2, bordering plot 1). Soil temperatures were measured at 5, 10, 20, 30 and 40 cm depth with copper-constantan thermocouples. Soil moisture was measured at 5, 20 and 40 cm depth (ML2x ThetaProbe, Delta-T Devices Ltd., Cambridge, UK). Precipitation was measured at 1 m height (rain and snow gauge with a heated funnel and base, Series 375, Met One Instruments, Inc., Grants Pass, OR, USA) at the center of the clear-cut. Four groundwater pipes were installed on each plot to allow for monitoring of groundwater level and for groundwater sampling. The pipes were drilled down to 2.5 m depth, 25 m North, East, South and West of the center of each plot ( Figure 2). Groundwater level was continuously measured (Druck PDCR 1830, GE-sensing, Billerica, MA, USA) at one location (L1) on plot 2 (25 m North of the center of the plot) from October 2010 onwards. The groundwater levels in all eight pipes were manually measured multiple times in 2010 and in 2012. From the DEM ( Figure 2) the differences in elevation at the ground surface between pipe L1 and the other seven pipes on plots 1 and 2 were calculated. This allowed for generating time series with daily average data on groundwater level for all eight pipes on plots 1 and 2 by assuming that groundwater levels at the other seven pipes were related to water table fluctuations at L1. These data were used to interpolate groundwater surfaces for plots 1 and 2 using kriging [73], and to estimate the distance between the groundwater level and the soil surface for each cell in the DEM (cell size 1 m 2 ). Daily average distances between the soil surface and the interpolated groundwater surfaces were calculated. The manual pipe measurements were used for validation of the estimated groundwater levels.
A short mast was installed centrally at the clear-cut site, close to the edge of plot 2, for measurements of net radiation (NR-lite net radiometer, Kipp and Zonen, Delft, The Netherlands) and air temperature. Air temperature was measured with a copper-constantan thermocouple in an aspirated radiation shield at 2 m height. Two web cameras (Mobotix M24M, Mobotix AG, Haupsitz, Germany) mounted on a 25 m tall mast took hourly digital photographs of both plots during daylight conditions (Figure 3), which allowed to track vegetation development and to detect when the ground was snow-covered. In addition, data on global radiation, incoming PAR, air temperature, air pressure, relative humidity and snow depth were further obtained from the Norunda main tower site, 600 m from the center of the clear-cut, either to complement existing measurements or to fill shorter gaps arising from power failures, etc. Large areas with standing water can be seen in the lower half of the image of plot 1. A small area with standing water can be seen in the lower left part of the image of plot 2, but this plot was in general drier than plot 1. The vegetation cover on both plots was dominated by grasses and shrubs, as most evident on plot 1.

Gradient Measurements
A 3 m tall tower was erected at the center of each plot and was equipped for measurements of turbulence parameters and concentrations of CO 2 , H 2 O, and CH 4 (May 2010-May 2013) as detailed below. Because of instrumental problems with the original gas analyzer, N 2 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. Data from the sonic anemometers were recorded on a computer. Air was continuously pulled (flow rates > 2.0 l min −1 ) from two heights at each tower through 100 m long high-density polyethylene tubes (inner diameter 8 mm) to a manifold inside a measurement hut placed centrally on the clear-cut. This setup where air was pulled from all intakes simultaneously to the manifold allowed for fast switching between the air intakes. The air intakes had 7 µm filters, were heated and had an orifice to create underpressure and prevent condensation of water vapor in the tubes. The intakes were placed at 0.6 m and 2.0 m height at the start of the measurements but the lower intakes were raised to 0.8 m in June 2012 to adjust for the increased height of the ground vegetation. The air stream from the intakes was analyzed for CH 4 , CO 2 and H 2 O with a Los Gatos tunable diode laser spectrometer (DLT-100, Los Gatos Research, San Jose, CA, USA) and N 2 O and H 2 O (QCL Mini Monitor, Aerodyne Research, Inc., Billerica, MA, USA) at a rate of 1 Hz. The data were recorded by a CR1000 datalogger (Campbell Scientific, Logan, UT, USA). Each location was sampled for 15 min each hour. These 15 mins were set up such that each level was sampled for 150 s before switching to the next level. The first 45 s of data from each level were discarded while the remaining data was used to calculate an average concentration for the two levels (N = 315 per level) during the 15-min period. Hence, an average concentration gradient for each of the gas species could be calculated once per hour and tower. These gradients were assumed to be representative for the full hour. Gas concentration data were quality controlled and filtered based on absolute limits wide enough to allow for high concentration values during stable stratification (300-800 ppm for CO 2 , 1.5-4.0 ppm for CH 4 , 0-30,000 ppm for H 2 O and 300-600 ppb for N 2 O). Gas concentration data were also filtered based on standard deviation following Schmid, Grimmond [74]: Any concentration value C i that deviated more than the product of a discrimination factor D and the standard deviation σ j from the mean concentration C j over the averaging interval was considered a spike and was removed. This was an iterative process where D increased from 3.5 to 5.0 in steps of 0.3.
Gas concentrations were corrected for dilution effects by water vapor. In addition, CH 4 and N 2 O concentrations were also corrected for dilution by CO 2 . All concentrations were converted from µmol mol −1 to µmol m −3 using the ideal gas law and measurements of air temperature and air pressure. A linear decrease in air pressure of 12.5 Pa m −1 with height was assumed in the conversions at the two measurement heights.

Eddy Covariance Data
The 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, version 5.0 (LI-COR, Lincoln, NE, USA), following standard eddy covariance methods. Turbulent fluctuations of the times series data were extracted using block-averaging over 15-min intervals corresponding to the gradient measurements. Raw data were de-spiked following Vickers and Mahrt [75] before calculating 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 [76], double rotation for tilt correction, high-pass filtering [77] and low-pass filtering [78]. Ogive tests [79,80] showed that 15-min averaging intervals were sufficient to capture most of the low-frequency contributions to the fluxes under turbulent conditions.

Flux-Gradient Calculations
Greenhouse gas exchanges between the ecosystem and the atmosphere were determined using an integrated eddy covariance and gradient method (flux-gradient (FG)). The method is based on Monin-Obukhov similarity theory and is described in detail by Simpson, Edwards [53] and Denmead [71]. Vertical turbulent fluxes F x of trace gases can be estimated as the product of the vertical concentration gradient of the gas species x times its turbulent diffusivity: where K x is the turbulent diffusivity (m 2 s −1 ) of gas x, C x is the concentration (µmol m −3 ) of gas x and z (m) is the measurement height. To account for non-neutral atmospheric conditions, K x can be expressed as where κ is the von Kármán constant (0.40, following Högström [81]), u * is the friction velocity, ϕ H is the universal function for sensible heat and ς is a dimensionless stability parameter. The universal function for sensible heat is commonly used for trace gases and latent heat [82] (p. 58). The stability parameter ς can be expressed as [82] (p. 55): where L is the Obukhov length, g is the gravity (9.81 m s −2 ), θ s is the characteristic potential temperature, ρ is the air density, C p is the specific heat capacity of air at constant pressure and H is the sensible heat flux. Equations (2) and (3) have to be integrated for an atmospheric layer between the air sampling heights z 1 (lower) and z 2 (upper), and adjusted for the zero-displacement height d. With that, F x can be expressed as where Ψ(ς) is the integrated form of the universal function [81]. The zero-displacement height d and roughness length z o were estimated from wind profile measurements during the summer of 2012. Five cup anemometers (cup array from ARI-49 (former USSR), mechanics and electronics developed by Tartu Observatory, Tartu, Estonia) were installed between 1 and 4 m height at a small mast placed centrally at the clear-cut. The cup anemometers were calibrated before and after the measurement period to make sure the obtained wind profiles were of high quality. The estimated d (0.2 m) and z o (0.07 m) were assumed to be valid for both plots throughout the measurement period (2010-2013) but were adjusted depending on the presence of snow (z o ) and on the thickness of the snow cover (d). The roughness length for a flat snow cover can be in the order of 0.001 m [83], but due to the small variations in topography, tree stumps, etc., at our plots, we chose a slightly higher value (0.003 m) for thick snow covers and an intermittent value between the two cases when the snow cover was thinner than the estimated d. The choice of d and z o does not affect the calculation of flux values but has an effect on footprint calculations (mostly z o ). See Table 1 for details. Figure 4a shows a land cover classification map of the clear-cut site, derived from a field survey of the plot treatment and from an airborne LiDAR survey in 2011. The two-dimensional footprint parametrization of Kljun, Calanca [84] was used for calculating hourly flux footprints for plots 1 and 2. From the footprint climatologies, derived by accumulating these hourly footprints for the whole measurement period, it can already be seen that the fluxes measured at each plot can be contributed to the corresponding land cover class to the largest part (up to 80%, see Figure 4b). Finally, merging hourly footprints with the land cover map allowed for estimating flux contributions per unit area (% per m 2 ) from the six land cover classes (class 1-4: plots 1-4; 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). See Figure 4 for details. Following Horst [85], the flux footprint modeling is based on the assumption that the flux footprint is representative for flux-gradient footprints obtained when the measurement height is defined as the arithmetic mean height of the two concentration measurements for stable stratification and is defined as the geometric mean height for unstable stratification.

Flux Data Filtering
Data quality tests on developed turbulent conditions and steady-state conditions were performed according to Mauder and Foken [86], resulting in a quality flag system (0-1-2) of the data. Sensible heat flux data with flag 2 were discarded from further analyses.
The errors of the gradient measurements were estimated to be 0.3 nmol mol −1 (CH 4 ), 0.15 µmol mol −1 (CO 2 ) and 4000 µmol mol −1 (H 2 O) for LGR data and 0.028 nmol mol −1 (N 2 O) and 1400 µmol mol −1 (H 2 O) for Aerodyne data, based on lab tests of the instruments RMS noise. Maximum errors in turbulent diffusivities were found to be 0.02 m s −1 for all conditions with u * > 0.06 m s −1 and hence, data collected during periods of u * < 0.06 m s −1 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 2 ), 20 µmol m −2 h −1 (CH 4 ) or 40 µg m −2 h −1 (N 2 O).
Careful examination of the flux data revealed a friction velocity dependency of the fluxes. The threshold where fluxes no longer increased with increasing u * varied slightly between seasons and between years but average u * thresholds of 0.16 ms −1 (plot 1) and 0.15 ms −1 (plot 2) were used for filtering out low-quality data during 2010-2013.
The GHG flux data were further filtered according to the flux footprints in two different ways. Firstly, hourly flux data were filtered according to the 80% footprint extent, ensuring a relatively homogeneous source area; i.e., fluxes with an 80% footprint extending outside the 50 m radius (black circles in Figure 4a) were discarded from further analysis. This approach minimized any possible contributions to the measured fluxes from other land-cover classes and allowed for a more direct comparison of the net GHG fluxes at plots 1 and 2, 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 classification, which, since it was based on the real plot geometry, maximized the plot areas (see Figure 4a 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 (all GHGs, 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.
By selecting data (filtered for footprint extent) measured during the same hour on both plots, possible problems with temporal autocorrelation of the fluxes were reduced and a non-parametric paired statistical test was carried out (Wilcoxon signed-rank test for zero median) in Matlab.

CO 2 Flux Data
Gaps in CO 2 flux data were filled using the R software (version 3.3.2) and the R package REddyProc [87]. REddyProc is a gap-filling and flux-partitioning tool that follows Reichstein, Falge [88]. 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 and Taylor [89] model is then used to derive an empirical model of R eco based on air temperature (T): where R eco,Tref is base respiration at reference temperature T ref (set to 10 • C), E 0 ( • C) is the temperature sensitivity, T 0 is kept constant at −46.02 • C as in Lloyd and Taylor [89]. E 0 and R eco,Tref are then estimated for each flux-averaging interval, resulting in hourly estimates of R eco . 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, Falge [88] and Wutzler, Lucas-Moffat [87] for a more detailed description of the gap-filling and flux-partitioning procedures.

CH 4 and N 2 O Flux Data
It was not possible to find strong relationships between methane exchange and different environmental variables at the clear-cut plots. Hence, data gaps 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 with outliers (σ > 5) removed prior to the averaging.
Gaps in N 2 O flux data were filled using linear interpolation of daily averaged data, similar to the method used by Scanlon and Kiely [90]. Daily averages were calculated in the same way as for the CH 4 data.

Global Warming Potential
For direct comparisons of the relative importance of the different GHGs, gap-filled CH 4 and N 2 O fluxes were also converted into CO 2 -equivalents by multiplication with a global warming potential factor in a 100-year perspective (GWP 100 ) of 28 (CH 4 ) and 265 (N 2 O) [91].

Vegetation Development
Landsat images with a spatial resolution of 30 m × 30 m were used to estimate vegetation cover on each of the experimental plots during the main growing seasons (June, July and August) of 2005 to 2014. The study site was captured by three stacks (paths 192-194; row 18) from Landsat 4, 5, and 7. Overlap among these three stacks was used to build up a high-resolution vegetation index time series for each of the experimental plots. High-quality data were identified by Fmask, which was used for automated detection of clouds, cloud shadows and snow masking for Landsat TM/ETM+ images [92,93]. The growing season Normalized Difference Vegetation Index (NDVI) [94] 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 m × 90 m areas overlapped the areas specified by the footprint classification ( Figure 4a) by 70-85%.

Environmental Conditions
The mean air temperature for the three study years (with a year defined as May 20 to May 19 of the following year in order to match the flux measurement periods) was 5.0 • C, 6.9 • C and 4. the first year than for the following years (610 mm and 611 mm, respectively). On average, soil moisture measured in organic soil, at 5 cm depth, was higher than in the mineral soil, at 20 cm depth. Figure 5 (panels d and e) suggests that a groundwater table within 20 cm from the soil surface did not result in saturated soil conditions at 20 cm depth in the soil moisture profile. This reflects the difficulties of representing a relatively large area with a highly disturbed soil surface, small-scale variability in topography, differences in soil properties, etc., by measurements at just one location. The groundwater table at the clear-cut 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 plot 2).
The soil surface at the groundwater pipe at the Norunda station was located at 45.5 m a.s.l., while the soil surface at the pipe on plot 2 was located at 49.3 m a.s.l. Thus, since the groundwater pipe in the forest was located 3.8 m lower than the pipe on the clear-cut, it is likely that the higher groundwater levels at the clear-cut were a consequence of the harvest. In general, clear-cutting can be expected to result in a raised groundwater level as a result of decreased canopy interception and evapotranspiration, but the exact magnitude and duration of this increase will depend on, for example, soil properties and vegetation development. During the second half of 2010, it became evident that plot 1 was wetter than plot 2, with several large areas with superficial groundwater. Plot 2 had some small areas with ponds ( Figure 3) but was, in general, much drier at the soil surface due to small differences in topography ( Figure 2) and possibly due to differences in soil properties (e.g., a clay layer has been found at approx. 1.5 m depth at some locations on both plots but the spatial extent of these layers is unknown). Standing water was often found in soil pits created during soil scarification (approx. 0.3-0.5 m deep) on both plots. The average distances between the calculated groundwater table and the soil surface for October 2010 to May 2013 on plots 1 and 2 ( Figure 6) were in good agreement with our qualitative view of the clear-cut, based on field experience, web camera images, etc. 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. The driest parts on plot 1 were located East of the tower, i.e., in the least represented wind directions. The wettest parts of plot 2 were in the SSE to NW wind directions, i.e., in the prevailing wind directions. On average, the groundwater table was closer to the surface on plot 1 than on plot 2 ( Figure 7). The daily average distance from the soil surface to the interpolated groundwater table on plots 1 and 2 were 0.16 m and 0.40 m, respectively, for the period October 2010 to May 2013. There was in general a good agreement between groundwater levels calculated for all pipes and the manual measurements (data not shown) when the groundwater level was close (>−0.2 m) to the surface (differences were less than 0.1 m). However, some of the extreme values (calculated groundwater table around 0.2 m above surface) following snowmelt in April 2011 and April 2013 ( Figure 7) are likely exaggerated, i.e., in reality, this surface water likely left the area as surface runoff.

Vegetation Development
The Landsat-derived NDVI for the clear-cut area ( Figure 8) suggests that there were no significant differences in vegetation cover on plots 1 and 2 before harvest in early 2009. After harvest, there was a sharp decline in summertime NDVI on both plots followed by a recovery phase until 2012-2013 where NDVI at both plots leveled off at a lower level than before harvest. In 2014, there were no longer any significant differences in NDVI between the plots. The recovery processes at plot 1 were faster and more substantial than on plot 2, indicating a higher vegetation cover from 2009 onwards. During 2009-2011, plot 1 had a significantly higher NDVI (p < 0.05) than plot 2. This is also supported by qualitative information from web cameras (see Figure 3). Given the low survival rate of the seedlings planted in 2010, most of the recovery in the summertime NDVI displayed in Figure 8 was due to the fast establishment of grasses and deciduous shrubs and trees.

Flux Footprints
The areas of plots 1 and 2 were, on average, well represented by the 80% annual flux footprints, despite hourly footprint variation with seasons and years (see Figure 4b for an example from 2012). The footprints were on average larger during nighttime, i.e., partly exceeded the plot area, and, as a result, more nighttime than daytime flux data were removed due to the footprint criteria.
Although the method of Horst [85] does not specifically require horizontally homogeneous fluxes, problems may arise if the concentration footprint of the air intakes at different heights in the towers are fundamentally different at a heavily disturbed clear-cut area. At a small spatial scale (<100 m 2 ), 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. However, at larger scales (>100 m 2 ), surface conditions were more homogenous because most disturbances at the surface were evenly distributed (e.g., mounding was done in a repetitive pattern, harvest residues ended up near the tree stumps). Nonetheless, e.g., wet patches with denser vegetation than on dry patches not seen to the same extent by both air intakes might have contributed to increased uncertainty in the measured fluxes and their source area.

CO 2 Fluxes
The average CO 2 fluxes (filtered for footprint extent; non-gap-filled) during the whole measurement period 2010-2013 were −0.3 µmol m −2 s −1 (plot 1) and 1.0 µmol m −2 s −1 (plot 2). The CO 2 fluxes at the two plots differed significantly (p < 0.001) for each individual year and for the whole period 2010-2013. On plot 1, the least represented wind directions were from the sector 30 • -120 • , an area containing some of the drier parts of this plot ( Figure 6). However, due to wet conditions and increased vegetation cover in the close vicinity of the tower, there was a net uptake of CO 2 (filtered for footprint extent) in the order of −1.2 to −1.5 µmol m −2 s −1 in this sector. For the remaining sectors, where footprints tended to include wetter areas, average fluxes (per 30 • sector) ranged from −0.2 to −1.5 µmol m −2 s −1 . On plot 2, the dominating wind direction was from 120 • -240 • (i.e., from the relatively wet parts of plot 2, see Figure 9a) but there was no clear dependence of wind direction on the CO 2 fluxes (filtered for footprint extent). Plot 2 was a net source of CO 2 with average net CO 2 emissions of 0.5-1.8 µmol m −2 s −1 in the different sectors.
Time series of gap-filled CO 2 data (daily averages, filtered according to footprint classification) are shown in Figure 9a. Initially, there was little vegetation on any of the plots and respiration dominated but already by the end of the first growing season, there were days with net uptake of CO 2 at plot 1. During the growing seasons of 2011 and 2012, the net uptake at plot 1 peaked in July, but by August, the fluxes were positive again. These early and short periods of net uptakes highlight the importance of grasses on clear-cuts for the CO 2 budget. The sudden shift in CO 2 fluxes on plot 1 in January 2013 was related to a raised groundwater table and increased soil moisture that possibly limited heterotrophic respiration and/or gross photosynthesis. Apart from a short period in June-July 2012, there were only daily average net emissions of CO 2 on plot 2, i.e., plot 2 was on a daily basis a source of CO 2 . Averaged over the whole period 2010-2013, the net emissions were 0.5 µmol m −2 s −1 (plot 1) and 1.1 µmol m −2 s −1 (plot 2), respectively. The cumulative CO 2 fluxes 2010-2013 were 3773 gCO 2 m −2 and 5181 gCO 2 m −2 on plots 1 and 2, respectively ( Table 2).  The difference in cumulative ecosystem respiration (R eco ) between the plots was relatively small (~600 gCO 2 m −2 for the three-year period), with the highest R eco on plot 1 (Figure 10). The lower (~1400 gCO 2 m −2 ) net emissions of CO 2 on plot 1 can be fully explained by higher GPP (~2000 gCO 2 m −2 ), which again highlights the importance of the fast establishment of grasses and shrubs on the total CO 2 budget. See Figure 10 for details on R eco and GPP during the individual years. Time series of daily average, gap-filled CH 4 fluxes (filtered for footprint classification) are shown in Figure 9b. Apart from a few periods with net CH 4 consumption, both plots were net sources of CH 4 . Plot 2 had relatively low and stable net emissions throughout the measurement period, while the CH 4 emissions at plot 1 were higher and more variable. During the first year of measurements, the net emissions of CH 4 were considerably higher on plot 1 than on plot 2, and there were several peaks (e.g., in August-September and in October) that coincided with precipitation events and peaks in soil moisture at 20 cm depth (  Figure 9c. The highest net emissions (up to~200 µgm −2 h −1 ) were on plot 1 during the summer of 2011, a relatively dry period with a few precipitation events causing large variations in groundwater level and soil moisture (especially at 20 cm depth). During December 2011 through mid-January 2012, wetter conditions were observed in general (at both 5 and 20 cm depth) and the groundwater table was often close to the surface. During this time, there was a small net uptake of N 2 O in the order of 5-10 µgm −2 h −1 on plot 1. The uptake period ended when the soil froze in mid-January 2012, and was instead turned into an emission peak that lasted for more than a month. During this period, soil moisture dropped rapidly when the soil froze and there was a decreasing trend in groundwater level. The N 2 O emissions started to decrease again in the middle of February when soil temperature was close to 0 • C and soil moisture levels and groundwater levels increased again. This was followed by a second emission peak at the beginning of March when the soil froze again. There was a second period with net N 2 O uptake during the latter half of April 2012, when soil temperature was above 0 • C again, the soil was wet and the groundwater level was close to the surface.
The net N 2 O emissions on plot 2 were generally higher than on plot 1, and also had a less variable temporal pattern. There were no periods with net uptake of N 2 O on the drier plot 2. The average N 2 O fluxes, based on the larger surface area specified by the footprint classification (Figure 4b (Figure 11, lower panel). The 2012 SON curve is based on less than a month of data and is thus not representative of the whole SON period, but is still interesting due to the scarcity of measurements of ecosystem N 2 O fluxes in managed forest ecosystems. There was a weak tendency towards larger diel variability in N 2 O exchange during periods with higher net N 2 O emissions and lower diel variations during periods with lower net emissions, or even net uptake (see Figure 9c). On the other hand, these low-emission periods were also characterized by emission peaks in

Total GHG Budgets
In order to compare the relative importance of the different GHGs, cumulative sums of gap-filled CH 4 and N 2 O fluxes were converted into grams of CO 2 -equivalents per square meter (in a 100-year-perspective). As shown in Table 2, the annual net CO 2 emissions were consistently lower on plot 1 than on plot 2, but this was partly offset by higher CH 4 emissions, especially during 2010-2011 and 2012-2013. During 2011-2012, N 2 O emissions contributed more to the total GHG budget on both plots than CH 4 emissions did, despite that N 2 O exchanges were only measured during 11 months this year.
The total net CO 2 and CH 4 emissions over three years were approx. 4107 gCO 2 -eq m −2 and 5274 gCO 2 -eq m −2 on plots 1 and 2, respectively. Carbon dioxide contributed 91.8% (plot 1) and 98.2% (plot 2) to the total CO 2 budget. Unfortunately, due to technical problems, the N 2 O fluxes were not measured for the full 36-month period, making it impossible to derive a full GHG budget for the whole period. In Table 2, we have added the contribution of N 2 O during the periods it was measured in the individual years. Despite the shorter measurement periods of N 2 O fluxes in the individual years, N 2 O still contributed significantly to the total GHG budget for these periods. To make it possible to directly compare the relative importance of the GHGs, we also compiled a total GHG budget for a full year, defined as 1 July 2011 to 30 June 2012 (see *** in Table 2), with simultaneous measurements of all GHGs. On plot 1, CH 4 and N 2 O contributed with 3.7% (39.3 gCO 2 -eq m −2 ) and 7.3% (77.9 gCO 2 -eq m −2 ), respectively, to the total GHG budget. On plot 2, CH 4 fluxes contributed with 1.5% (25.8 gCO 2 -eq m −2 ) while N 2 O fluxes contributed with 7.6% (128.8 gCO 2 -eq m −2 ) to the GHG budget.

Discussion
The results show that clear-cutting of a hemiboreal forest stand resulted in large GHG emissions over the first three years and that the relative contribution of the three GHGs in the total GHG budget depended on the degree of wetness of the plots. The clear-cutting resulted in waterlogged soils, which most likely caused the site to switch from a sink to a source of CH 4 . N 2 O fluxes contributed more to the total GHG budget than CH 4 fluxes did, even on the wet plot 1 where occasional net uptake of N 2 O was found.

Vegetation Development
The vegetation development had a large impact on CO 2 fluxes at both plots. While the ecosystem respiration was similar at both plots, the higher GPP (due to higher vegetation coverage) at plot 1 resulted in lower total net emissions than at the drier plot 2. This observation was also supported by both the summertime NDVI data (Figure 8) and by the differences in the temporal evolution of the diel curves in CO 2 fluxes ( Figure 11) on plots 1 and 2. Given the low survival rate of the seedlings planted in May 2010, it was evident that a large part of the differences in GPP between the plots could be explained by the faster and more substantial establishment of grasses and shrubs at the wetter plot 1. Even though GPP increased 5.4-fold from year 1 to year 3 on plot 2, it was consistently lower than at plot 1, where GPP only increased by a factor of~2.

Gap-Filling of CH 4 and N 2 O Fluxes
In order to derive a full GHG balance, gap-free data of high quality are required. Well-established gap-filling methods for times series of CO 2 fluxes exist; however, there is currently no general consensus on gap-filling of CH 4 and N 2 O fluxes. Attempts have been made to fill gaps in CH 4 fluxes of wetlands using, e.g., temperature relationships, simple linear interpolation and neural networks [95][96][97]. Sundqvist, Persson [98] described an empirical model based on multiple stepwise regression analysis for upscaling methane exchange from point measurements (chambers) to the stand scale at the Norunda forest using relationships between CH 4 fluxes and soil temperature, soil moisture and water table depth. It was not possible to find any similar relationships between CH 4 exchange and typical controlling variables on the clear-cut plots 1 and 2, most likely due to the heavily disturbed soil and changes in hydrology following harvest, and hence, we filled gaps using linear interpolation of daily average CH 4 fluxes. N 2 O fluxes are notoriously variable, both temporally and spatially, making gap-filling even more complex than for CH 4 . Scanlon and Kiely [90] used linear interpolation of daily average fluxes to fill gaps in N 2 O fluxes (EC) at a managed grassland, while Mishurov and Kiely [99] tested linear extrapolation, moving average and a look-up table approach to fill gaps of different lengths (day, month and year) in the same ecosystem. In total, six methods were tested by the authors and they were found to agree reasonably well. Additionally, Merbold, Eugster [31] filled gaps in N 2 O fluxes at a managed grassland using look-up tables for periods with similar environmental conditions. For our highly disturbed ecosystem with large temporal and spatial variability in known controlling parameters, and poor spatial representation of, e.g., soil moisture and soil temperature measurements, the simpler linear interpolation method was chosen.

Diel Patterns of CH 4 and N 2 O Fluxes
The absence of clear diel variations in CH 4 exchange at the drier plot 2 was probably a consequence of the, on average, larger distance from the soil surface to the groundwater, which resulted in lower CH 4 production rates. In addition, the highly disturbed soil surface is likely to have reduced CH 4 consumption, as compared to undisturbed soil. Teepe, Brumme [100] and Strömgren, Hedwall [28] showed decreased CH 4 consumption, and even net CH 4 emissions, from wheel ruts on clear-cuts, probably as a consequence of reduced macropore volume due to compaction, which limited diffusion of atmospheric CH 4 into the soils. Other studies have shown that it takes years [101] to more than a century [9] for disturbed soils to regain their full sink capacity.
There were either clear or weak diel variations in CH 4 exchange on the wetter plot 1 during all seasons, but the amplitude varied between seasons and years. Given that methanogens normally have a higher temperature sensitivity than methanotrophs [102], it is possible that the high groundwater level and the severely disturbed soil surface (i.e., high production rates and low consumption rates) can explain the observed diel patterns with daytime maxima in CH 4 exchange. On the other hand, we found no, or only very weak, correlations between soil temperature (at any depth) and CH 4 fluxes at any of the plots. However, this might be a result of the poor spatial representation of soil temperature measurements at the clear-cut. The absence of a diel pattern on plot 1 during JJA in 2011 was probably related to the lower than average groundwater level during this period.
A few studies have reported diel patterns of CH 4 exchange in forest ecosystems. Pattey, Strachan [48] reported daytime maxima of CH 4 emissions over a boreal forest, using both flux-gradient and EC techniques during shorter campaigns, and noted that the CH 4 cycle closely followed the soil temperature cycle at 5 cm depth. Other studies have reported daytime minima [49,103] in CH 4 fluxes. Sundqvist, Mölder [104] measured CH 4 fluxes in the Norunda forest with different gradient methods above the canopy and suggested that the observed daytime minima could be a result of plant uptake (as previously observed from branch chamber measurements by Sundqvist, Crill [105]) and that the observed nighttime maxima could be a result of larger nighttime footprints, with distant, wetter source areas (such as the clear-cut in this study) contributing more to the measured fluxes than during daytime.
Diel patterns of N 2 O fluxes from agricultural soils and grasslands have been known for about 40 years [106][107][108]. Recently, Shurpali, Rannik [109] reported clear diurnal variability in N 2 O fluxes (measured with EC on agricultural land) depending on soil N availability. They attributed daytime flux maxima during high soil N availability (i.e., following fertilization) to soil temperature control on the involved microbial processes, while the N 2 O fluxes during low soil N availability correlated negatively with GPP. Their study suggests that plants play a role in the N 2 O exchange by moderating soil O 2 levels in the soil and nitrogen availability (either directly through plant uptake of nitrogen or indirectly through root exudates of carbon that stimulates immobilization of inorganic nitrogen by heterotrophic bacteria). The net result would be reduced N 2 O production and increased reduction of N 2 O to N 2 and hence, reduced net emissions or even net uptake of N 2 O. Our results are different from those of Shurpali, Rannik [109], i.e., with clear daytime maxima of the N 2 O exchange during seasons when GPP was at its maximum, and during what would be considered as a low-flux period in the study of Shurpali, Rannik [109]. During seasons with no or low diel amplitude of the diel CO 2 exchange (DJF, MAM), the amplitude of the diel N 2 O exchange was also non-existent or low. During seasons with larger diel amplitudes in CO 2 fluxes (JJA, SON), the diel amplitude in N 2 O fluxes was also higher. Furthermore, the magnitude of the flux seemed to have an impact on the amplitude of the diel curves; the lower the flux, the lower the amplitude. Zona, Janssens [110] measured EC N 2 O fluxes above a short-rotation poplar plantation and found no diurnal patterns in the fluxes except during shorter periods. The controlling factor of the diurnal pattern during these short periods was found to be wind speed.
Unfortunately, we do not have any data on soil O 2 concentrations or on available soil N that could explain our observed diel patterns at plots 1 and 2. However, it is likely that with flux footprints integrating fluxes from a heavily disturbed ecosystem with small-scale variability in topography, soil moisture, soil temperature, distance to the groundwater table, soil properties, pH, available soil N, etc., hence encompassing a wide range of environmental, biogeochemical and biogeophysical conditions, it is difficult, if not impossible, to draw any conclusions on the key drivers of the observed diel patterns.
We are not aware of any micrometeorological studies of N 2 O fluxes at clear-cut forest stands.

Greenhouse Gas Fluxes
Our approach, applying a very strict filtering based on footprint extent (and on turbulence parameters, etc.) of the flux data resulted in a much reduced dataset and in a strong confidence that the remaining fluxes were representative for plots 1 and 2, respectively. By applying a statistical test (Wilcoxon signed-rank test for zero median) on pairwise data (fluxes measured during the same hour) from plots 1 and 2, robust and statistically significant differences in fluxes of CO 2 , CH 4 and N 2 O were observed. We mainly attribute the differences in fluxes to differences in the wetness at the plots, either directly (e.g., through the effect of the groundwater table position on fluxes of CH 4 and N 2 O) or indirectly, through the faster and more substantial establishment of vegetation on the wetter plot 1 (and especially around ponds and other wetter areas at plot 1), which had a significant impact on the net emissions of CO 2 . As a consequence of the very strict data filtering on footprint extent, a large proportion of nighttime and wintertime data was lost and hence, the non-gap-filled CO 2 fluxes were biased towards net uptake of CO 2 (plot 1) or smaller net emissions (plot 2).
The less strict filtering on footprint classification, and the subsequent gap-filling, made it possible to derive annual cumulative sums of NEE, GPP and R eco for plots 1 and 2. The values for the first three years after clear-cutting reported here (Table 2) were lower than those reported after clear-cutting of Douglas fir stands in Western Canada Humphreys, Andrew Black [40] and Paul-Limoges, Black [45] reported total net emissions of CO 2 of~6380 gCO 2 m −2 and~8900 gCO 2 m −2 , respectively). Additionally, Korkiakoski, Tuovinen [15] reported higher net emissions from a drained peatland in Finland for the first two years following clear-cutting of a Scots pine-dominated stand (~3090 gCO 2 m −2 and 2070 gCO 2 m −2 for year 1 and 2, respectively). Our net CO 2 emissions reported here are higher than those reported by Coursolle, Giasson [43] (net emissions of~1550 gCO 2 m −2 the first three years after clear-cutting of a Black spruce (Picea mariana) in eastern Canada). Amiro, Barr [33] reported annual net emissions of CO 2 in the order of 730-3120 gCO 2 m −2 following clear-cutting of a wide range of forest across North America, with boreal forests being in the lower end of that range.
Similar to Amiro, Barr [33] and Williams, Vanderhoof [44], we observed relatively stable R eco for the years following clear-cutting and a gradual recovery in GPP, which further stresses the importance of vegetation development for the carbon balance of a clear-cut forest stand.
We do not have any pre-harvest measurements of CH 4 fluxes at our clear-cut but simultaneous chamber measurements in an undisturbed forest in the close vicinity of the Norunda station and on plot 2 at the clear-cut by Sundqvist, Vestin [14] showed net uptake in the order of −10 µmol m −1 h −1 in the undisturbed forest and average net CH 4 emissions of 13.6 µmol m −1 h −1 on plot 2 in October-November 2010. The results from plot 2 were in good agreement with our flux-gradient CH 4 measurements during the same period. Sundqvist, Vestin [14] used an automated chamber system with five chambers placed on both disturbed and undisturbed soil on plot 2 and found net emissions in the order of 15-32.5 µmol m −1 h −1 from chambers placed on disturbed soil and either small net emissions or small net uptake (−3 µmol m −1 h −1 ) from chambers placed on undisturbed soil. The distance from the chamber frame to the groundwater table was 0.44-0.5 m for the chamber with net uptake, while the distances were less than 0.21 m for chambers with net emissions. Most chambers were placed southwest of the tower on plot 2, i.e., at the wetter parts of plot 2 ( Figure 6). The average distances to the groundwater table were 0.16 m and 0.40 m on plots 1 and 2, respectively, which most likely explains the observed differences in CH 4 emissions between the plots. Furthermore, our observations during the period 2010-2013 support the initial observation by Sundqvist, Vestin [14] in 2010 that waterlogging following clear-cutting caused the site to switch from a sink to a source of CH 4 . Others have also reported how a site switched from a sink to a source following clear-cutting, either as a result of increased soil moisture that inhibited CH 4 consumption [13], a higher water table [15] or as the combined result of increased soil moisture, higher water table and increased soil N availability [16].
Our results with high and variable CH 4 emissions (plot 1) or low and steady emissions (plot 2) are different from results reported by Strömgren, Hedwall [28], who used chambers to measure fluxes from soils disturbed by soil preparation and wheel ruts in three clear-cut Norway spruce forests in central Sweden. The authors found reduced uptake rates of CH 4 at soil pits and large emissions from wheel ruts, as compared to control plots with undisturbed soils. However, when scaled to stand level according to the area coverage by each disturbance type, they found no significant differences between control plots and treated plots. Similarly, Kulmala, Aaltonen [68] found no significant difference in CH 4 fluxes between a mature forest and a clear-cut plot in a chamber study at a Norway spruce forest in southern Finland. The soil at the clear-cut remained a sink throughout the three-year study period.
The net fluxes of N 2 O at plot 1 were highly variable, ranging from a small net uptake (December 2011-January 2012 and April-May 2012) to large emission peaks (e.g., JJA 2011, February 2012 and March-April 2012). It seems likely that the emission peaks in JJA 2011 were caused by precipitation events and related fluctuations in soil moisture and groundwater table depth and that the emission peaks in February 2012 and March-April 2012 was related to freezing and thawing of the soil at plot 1. N 2 O fluxes are known to be highly variable in space and time and many have reported that a single, or a few, emission events following, e.g., fertilization [111], precipitation events [110], or freeze-thaw events [112], can make up a large portion of the annual N 2 O budget. The N 2 O fluxes at plot 2 were generally higher and less variable than at plot 1. There were no clear effects of precipitation events and freeze-thaw cycles at plot 2, possibly because of the lower groundwater table.
The small net uptake of N 2 O we observed at the wet plot 1 is similar to previously reported EC values for forests [24,25,53] and for chamber measurements in the forest [112] and in clear-cut boreal forests [27,28].

GHG Budgets
For the only full year with simultaneous measurements of all GHGs (1 July 2011 to 30 June 2012), CO 2 contributed 89.0% and 90.9% to the total GHG budget on plots 1 and 2, respectively. The contribution of CH 4 to the total budget was low during this year, possibly as a consequence of the lower than average groundwater table position (during July and August in 2011 and during January and February 2012) and because of soil frost during the snow-free winter 2011-2012. The contribution of N 2 O to total GHG emissions was 7.3% and 7.6% on plots 1 and 2, respectively, although the N 2 O emissions at the drier plot 2 were 65% higher in absolute terms than at the wetter plot 1.
The total GHG budget was lower (plot 1) and higher (plot 2) than values reported (1581 gCO 2 -eq.m −2 with GWP 100 factors of 23 and 296 for CH 4 and N 2 O, respectively) for a chamber study during the first ten months following clear-cutting of a Sitka spruce plantation in Northern England [16]. Zona, Janssens [54] reported a net uptake of CO 2 of −76 gCO 2 -eq.m −2 and a combined net release of 351 gCO 2 -eq.m −2 for CH 4 and N 2 O (using GWP 100 factors of 25 and 298 for CH 4 and N 2 O, respectively) during a one and a half year study period at a short-rotation bioenergy poplar plantation in Belgium. In their study, water availability was the main driver of CO 2 uptake and CH 4 and N 2 O emissions. Peichl, Arain [32] combined EC measurements of CO 2 (mean value over a four year study period) with chamber measurements of CH 4 and N 2 O (mean value over two years) in a White pine (Pinus strobus) stand that was seven years old at the start of the measurements. As in our study, CO 2 fluxes were the main contributor to the total GHG budget, while CH 4 and N 2 O fluxes combined contributed 13%. Korkiakoski, Tuovinen [15] also used a combination of EC measurements of CO 2 with manual chamber measurements of CH 4 and N 2 O at a clear-cut forest on a nutrient-rich drained peatland and noted a switch from a weak CH 4 sink to a weak CH 4 source and from N 2 O neutral into a significant N 2 O source. While the CO 2 emissions were large, the CH 4 emissions were deemed negligible in a global warming potential perspective (GWP 100 = 34) and the N 2 O contribution was in the order of 10% (GWP 100 = 298) of the CO 2 fluxes, i.e., very similar to our study.
At plot 2, CH 4 emissions were relatively low and constant during the study period (Table 2), while the net emissions were much lower at plot 1 during 2011-2012, which to a large part overlaps with the year used for creating full GHG budgets (see Table 2). During the periods 2010-2011 and 2012-2013, CH 4 fluxes (in gCO 2 -eq.m −2 ) amounted to 9.5% and 12.9% of the CO 2 emissions at the wetter plot 1 and only 3.6% during 2011-2012. At the drier plot 2, the CH 4 emissions never amounted to more than 2.3% of the CO 2 emissions.
This study suggests that forest management strategies aiming at enhancing the carbon sink strength in forests also need to consider the impact of forest management on the fluxes of CH 4 and N 2 O since these can contribute with approx. 10% each of the measured CO 2 emissions during individual years. From the 20-year perspective, the importance of CH 4 fluxes (GWP 20 = 84) increases three-fold, which has to be taken into account for any management strategy aiming at maximizing the mitigation potential of forests in the near future.
In addition, the impact of a future, possibly more variable climate with increased frequency and intensity of precipitation events, more extreme temperatures and droughts on GHG budgets should also be considered when implementing forest management policies and practices aiming at enhancing the mitigation potential of forests.

Conclusions
We measured the fluxes of CO 2 , CH 4 , and N 2 O using a flux-gradient technique at a clear-cut hemiboreal forest stand during three consecutive years, starting one year after the clear-cutting. The method proved to be a reasonable compromise between cost (of instruments, maintenance, running costs, etc.), data coverage and accuracy, and was well suited for inter-plot comparisons even in a heavily disturbed ecosystem. However, future studies should address the large spatial variability in soil properties, soil temperature, soil moisture, etc., as well as soil biogeochemistry to allow more rigorous and robust gap-filling strategies for CH 4 and N 2 O fluxes and to make it possible to obtain a better insight into which processes control the GHG fluxes at clear-cut forest sites. Ideally, future studies should also include measurements to disentangle the relative importance of production and consumption of GHGs, and their drivers, in soils at clear-cut forest stands in order to facilitate modeling and upscaling efforts.
We observed significant emissions of GHGs both at a wet and a dry plot. The relative contribution of CO 2 , CH 4 and N 2 O to the GHG budget differed depending on the degree of wetness, which affected GHG fluxes both directly and indirectly. The faster establishment of vegetation on the wet plot resulted in lower net emissions of CO 2 but this was partially offset by significantly higher CH 4 emissions than at the dry plot 2. It seems likely that waterlogging after clear-cutting caused both plots to switch from a net sink to a net source of CH 4 . We observed significant N 2 O emissions at both plots. The short periods with a low net N 2 O uptake at the wet plot 1 were likely a consequence of an increased reduction of N 2 O to N 2 by denitrifying bacteria during periods when the groundwater table was at, or close to, the surface.
Any policy aiming at enhancing the mitigation potential of forests should consider the combined effect on the climate system of GHG fluxes following forest management operations and also consider what impact a future, more variable climate might have on GHG budgets of forests.