Effects of clearfell harvesting on soil CO2, CH4 and N2O fluxes in an upland Sitka spruce stand in England

The effect of clearfell harvesting on soil greenhouse gas (GHG) fluxes of carbon dioxide (CO2), methane (CH4) and nitrous oxide (N2O) was assessed in a Sitka spruce forest growing on a peaty gley organo-mineral soil in northern England. Fluxes from the soil and litter layer were measured monthly by the closed chamber method and gas chromatography over four 10 years in two mature stands, with one area harvested after the first year. Concurrent measurements of soil temperature and moisture helped to elucidate reasons for the changes in fluxes. In the three years after felling, there was a significant increase in the soil temperature, particularly between June and November (3 to 5 C higher), and in soil moisture which was 62% higher in the felled area, and these had pronounced effects on the GHG balance, in addition to the removal of the trees and their carbon input to the soil. Annual soil CO2 effluxes reduced to almost a third in the first year after felling (a drop from 24.0 to 15 8.9 t CO2 ha yr) and half in the second and third year (mean 11.8 t CO2 ha yr) compared to before felling, while those from the unfelled area were little changed. Annual effluxes of N2O more than doubled in the first two years (from 1.0 to 2.3 and 2.5 t CO2e ha yr, respectively), although by the third year they were only 20% higher (1.2 t CO2e ha yr). CH4 fluxes changed from a small net uptake of -0.03 t CO2e ha yr before felling to a small efflux increasing over the 3 years to 0.34 t CO2e ha yr, presumably because of the wetter soil after felling. 20 Soil CO2 effluxes dominated the total annual net GHG emission when the three gases calculated were compared using their global warming potential (GWP) of the three gases, but N2O contributed up to 20% of this. of the total GWP annual emissions. This study showed fluxes of CO2, CH4 and N2O responded differently to clearfelling due to the significant changes in soil biotic and abiotic factors and showed large variations between years. This demonstrates the need for multi-year measurements of all GHGs to enable a robust estimate of the effect of the clearfell phase on the GHG balance of managed 25 forests. This is one of a very few multi-year monitoring studies to assess the effect of clearfell harvesting on soil GHG fluxes.

the absence of evapotranspiration from trees (Zerva and Mencuccini, 2005;Wu et al., 2011;Kulmala et al., 2014;Sundqvist et al., 2014;Korkiakoski et al., 2019); soil temperature through reduced shading (Zerva and Mencuccini, 2005;Wu et al., 2011;Kulmala et al., 2014); soil bulk density, due to soil disturbance and compaction caused particularly by mechanised harvesting equipment (Yoshiro et al., 2008;Mojeremane et al., 2012). Chemical factors include: soil pH (Smolander et al., 1998;Kim, 2008;Kulmala et al., 2014;Sundqvist et al., 2014) and inputs of soil N and organic C (Smolander et al., 2015;45 Tate et al., 2006;Hyvönen et al., 2012;Kulmala et al., 2014;Sundqvist et al., 2014), due to nutrient release from the decomposition of residual organic matter and root biomass. Interactions between these factors, particularly with the loss of plant litter input and root activity after felling, will strongly affect soil biological processes responsible for production and consumption of GHGs. They will particularly affect nitrous oxide (N2O) production from aerobic nitrification and anaerobic denitrification, methane (CH4) production by methanogenic organisms from anaerobic decomposition in oxygen-poor 50 environments or uptake by methanotrophs through oxidation in aerated soils, and carbon dioxide (CO2) efflux during respiration and decomposition and uptake during photosynthesis. The cessation of the autotrophic respiration (Ra) component of the total soil respiration (Rt) after felling should cause a large decline in CO2 efflux as a meta-analysis of soil respiration partitioning studies reported that the Ra/Rt ratio in temperate forest soils ranges from 20 to 59% (Subke et al., 2006). In addition the death of tree roots after felling will inhibit the microbial decomposition of root exudates reducing the CO2 effluxes further. 55 Conversely, the increased soil temperatures following tree removal and higher nutrient availability from decomposition of litter and other plant materials such as brash may increase soil heterotrophic respiration (Rh) (Yashiro et al., 2008).
There are a wide range of results published on the effects of forest clearfelling on soil GHG fluxes. For CO2, the overall soil efflux after felling depends on the balance between the decrease in root respiration and the increase in microbial respiration in response to warmer soil temperatures and nutrient availability from decomposition of litter and other plant 60 materials such as brash. Therefore, sSome studies have shown CO2 effluxes were reduced (Striegl and Wickland, 1998;Laporte et al., 2003;Zerva and Mencuccini, 2005;Goutal et al., 2012, Korkiakoski, 2019, increased (Tate et al., 2006;Kim, 2008;Kulmala et al., 2014), or unchanged (Toland and Zak, 1994;Butnor et al., 2006;Takakai et al., 2008;Yoshiro et al., 2008).
The differences between studies on the effect of harvesting disturbance on GHG fluxes have also been attributed to microsite distribution within the harvested area (Laporte et al., 2003). According to Lavoie et al. (2013), the impact is often site specific, 65 affected by the severity of the disturbance or removal of surface organic matter (Fleming et al., 2006;Moroni et al., 2007;Slesak et al., 2010) or by the length of time following harvest (Humphreys et al., 2006;Peng and Thomas, 2006;Olayuyigbe et al., 2012). Goutal et al. (2012) examined the duration of physical, chemical and biological disturbances in the soil following mechanized harvesting of a beech/oak forest in NE France. Their measurements showed that soil CO2 effluxes reduced which they attributed to an increase in the frequency and duration of anoxic conditions resulting from poor soil gas diffusion after 70 heavy forestry traffic. Similarly, Kulmala et al. (2014) observed a slight decrease in soil CO2 efflux in the first growing season after clear-cutting of a boreal Norway spruce stand, although this was probably due to decreased tree root respiration. However, during the following two years, CO2 efflux at their clear-cut site was significantly higher than in their mature stand, which was attributed to increased decomposition, stimulated by higher soil moisture and temperature. They observed no significant difference in CH4 uptake due to clear-cutting. However, others have shown that after harvesting forest soils turned from a CH4 75 sink to a source (Zerva and Mencuccini, 2005;Castro et al., 2000;Takakai et al., 2008;Yashiro et al., 2008;sSundqvist et al., 3 Therefore, there is an urgent need to understand and quantify with long-term measurements the effect of forest clearfelling on the GHG budget and to incorporate these into life cycle analyses for the complete forest growth cycle for the main forest 85 systems relevant to a country or region (Skiba et al., 2012). Therefore,Tthe objective of this study was to conduct a relatively long-term (4 year) assessment of the effect of clearfell harvesting on soil (including litter layer) fluxes of CO2, CH4 and N2O in a spruce plantation on organo-mineral soil that is typical of many British upland forests.

Site description and study layout 90
The study site was in Harwood Forest, (55°13.1' N, 2° 1' W), Northumberland, north-east England, and comprised a forest area of approximately 4000 ha with an elevation of 200 to 400 m above sea level (Fig. 1). The regional climate is temperate oceanic, with a mean annual rainfall of 1472 mm and mean air temperature of 7.5°C (min. -7 and max. 26.4°C), measured by an automatic weather station (AWS) mounted above the tree stand during the period from April 2015 to April 2018. The tree cover consisted predominantly of even-aged Sitka spruce (Picea sitchensis, Bong. Carr.) stands. Sitka spruce is the most 95 common conifer in Great Britain (GB) and represents 26% of the total forested area in GB and 51% of the conifer area (Forestry CommissionForest Research, 2020. The main soil type is a seasonally waterlogged organic-rich peaty gley classified as histogleysol with a peat (O Horizon) thickness varying from 15 to 4s0 cm (Zerva and Mencuccini, 2005), developed in clayey glacial till derived from carboniferous sediments (Pyatt, 1970). For the purpose of this study, two nearby areas (A and B, 2 km apart) of mature stands within the forest with similar previous management, soil type and elevation (280 to 290 m) were chosen to carry out measurements over four years between 18 February 2014 to 16 April 2018. During this period, one area (A) was left unfelled (A-yr1 to A-yr4), and the other (B) was 110 clearfelled after one year (B-yr1 before felling and B-yr2 to B-yr4 after felling). In area A, the 40-ha stand was of second rotation, even-aged mature Sitka spruce planted in 1973, with yield class of 18 m 3 ha -1 and mean tree density of 1348 trees ha -1 . In area B the Sitka spruce stand was planted in 1958, with yield class of 16 m 3 ha -1 and mean tree density of 1375 trees ha -1 prior to felling, and the felled area covered 42 ha. Felling operations were carried out between late January and early March 2015 and followed standard practices of the Forestry Commission (Murgatroyd and Saunders, 2005). Only timber larger than 115 7 cm diameter was removed from site, leaving tree tops and branches on site in rows, with some used as brash-mats to prevent compaction of soil by the heavy harvesting machinery.

Gas flux measurements and analysis
Forest soil fluxes of CO2, CH4 and N2O were measured at approximately monthly intervals over the four-year study period from the two areas of the forest. Flux measurements were made using a modified design of the manual static chamber method 120 described by Yamulki et al. (2013). Each chamber was constructed of opaque PVC with dimensions of 40 cm × 40 cm × 25 cm height to provide a volume of 40 L and placed temporarily on permanently installed frames. The frames were inserted tightly into the ground to a depth of about 3 cm prior to the start of the measurements. The bottom of the chamber had a neoprene rubber foam gasket to ensure a gas-tight seal with the frame, and the top of the chamber had a pressure vent.
Within each area, 8 chambers were positioned randomly in a transect within a 100 m 2 area. During each gas flux 125 measurement, the chambers were placed on top of the frames for up to 60 min and duplicate gas samples of the chamber headspace were taken immediately after closure and then at 3 subsequent 20-min intervals. Gas samples were taken after the chamber was closed by connecting a polypropylene syringe to a chamber sampling port fitted with a three-way stopcock. The syringes were immediately used to fill (under atmospheric pressure) pre-evacuated 20 mL vials fitted with chlorobutyl rubber septa. Concentrations of CO2, CH4 and N2O were determined within a week using a headspace-sampler (TurboMatrix 110) 130 and gas chromatograph (GC, Clarus 500, PerkinElmer) equipped with an electron capture detector (ECD) for N2O analysis, a flame ionization detector (FID) for CH4 analysis, and a catalytic reactor (methanizer) for CO2 analysis by reducing CO2 to CH4 before analysis by the FID detector. The repeatability of the GC gas analysis (assessed as 3 × standard deviation of 20 repeated measurements of standard CO2, CH4 and N2O concentrations at ambient levels) was better than 4% for all gases.
Gas fluxes were calculated based on linear increases of gas concentrations inside the chambers with time. For CO2 135 however, if the concentration increase was not linear then fluxes were determined using the R HMR package (version 1.0.1) to plot a best-fit line to the data (Pedersen et al., 20210;R Core Team, 2016) to correct for the non-linearity. We did not apply the HMR model for CH4 and N2O fluxes as the non-linear fitting plot is very sensitive to variability and outliers in the measured GHG concentrations, particularly for low fluxes and with only four data points per chamber (as also noted by Pihlatie et al., 2013;Brümmer et al., 2017;Korkiakoski et al., 2017) which results in large apparent 'spikes', failure to calculate the fluxes 140 on many occasions, and likely overestimation of calculated GHG fluxes (Pavelka et al., 2018). If CO2 concentration changes with time were not significant, fluxes for all gases were rejected for that sample as this was judged to be indicative of gas leakage within the chamber headspace. Overall 18 samples were rejected, 9 of which were during snow fall period in January

2016.
Although most studies now measure soil CO2 fluxes with IRGAs (as noted by Yashiro et al., 2008), which is viewed 145 as more accurate than gas sampling and GC analysis, we needed to use the GC method to measure all 3 GHG within the logistical constraints of the experiment. Therefore, to give further assurance in the gas flux calculations, CO2 effluxes were also compared with those from another 25 static chambers (20 cm diameter, 4.2 L volume, LI-8100-103 Survey Chamber, Li-Cor Inc., Lincoln, Nebraska, USA) measured in situ in both areas by a closed loop infra-red gas analyser system (IRGA, LI8100A, Li-Cor Inc.), each over a 2-minute duration after chamber closure (Xenakis et al., 2021). These chambers were 150 positioned over a much wider area than the GC chambers to characterise spatial variation and effluxes were measured for three years but only after felling from Feb 2015. The results (Fig. 2S1) showed a mean flux difference over the 3-year period of only 6.8% higher by the IRGA method in the unfelled area and 19.5% higher in the felled area compared with the fluxes measured by GC. These are relatively small differences between the two methods, when considering the higher site heterogeneity in the felled area, the inherent differences in the analytical methodologies and as the vegetation in the IRGA chambers was not cut 155 as it was in the GC chambers. Figure. 2. Comparison between mean CO2 effluxes measured by gas chromatograph (GC) from the 8 static chambers in each forest area (A = mature stand, B = clearfell) and mean CO2 effluxes measured by infra-red gas analyser (IRGA), from 25 static 160 chambers in the same areas.
There was no understorey vegetation in the forest stands and no visual evidence of vegetation growth within the chambers in the first two years of the study, with the chamber ground surfaces covered by dense Sitka spruce needle litter.
However, in the third year there was growth of Juncus sp. in one chamber in the felled area B between May and September 165 2016 (maximum height was 40 cm before cutting but the total volume was < 5% of the chamber volume) so the Juncus was cut continually thereafter. The effect of the vegetation cut on the fluxes was assessed from the flux measurements in all the 8 chambers before and directly after first cutting. The statistical analysis revealed no significant differences between mean chamber fluxes before and after cutting for all gases indicating that variations between the chamber fluxes were greater than that due to cutting. In the third and fourth years there was a proliferation of moss in two of the chambers in area A. This could 170 not be removed without substantial disturbance of the soil surface; comparison of chambers with and without moss showed no significant differences between mean fluxes. For CO2, the effluxes measured will therefore be from aerobic and anaerobic decomposition processes, and respiration of soil organisms and roots.

Soil moisture and temperature measurements
During each flux sampling day soil temperatures ( o C) at 2 cm and 10 cm depths were recorded from one point around each 175 chamber and volumetric moisture content (m 3 m -3 ) at 6 cm depth was recorded from three points around each chamber. Soil temperature was recorded by a digital temperature probe (Hanna model Checktemp 1) and the volumetric moisture content by a moisture sensor (SM 200 attached to a handheld HH2 moisture meter, Delta-T Devices Ltd, Cambridge, UK). Sensitivity of CO2 efflux to temperature was determined with a Q10 function (Atkin et al., 2000 Shi andJin, 2016) which is the proportional change in respiration resulting from a 10°C increase in temperature, derived from data for each year using mean daily values 180 of CO2 efflux for each year using the equation: where b is the slope of the exponential relationship between soil CO2 effluxes and soil temperature at 2 cm depth, estimated 185 from a log-linear fit. The apparent Q10 values were calculated from the temperature measured near the soil surface (at 2 cm) as recommended by Pavelka et al. (2007), because the strength of the CO2 efflux and temperature relationship decreases with soil depth.

Soil parameters and root biomass measurements
To assess changes in soil characteristics due to the felling, additional soil parameters that were likely to have an impact on 190 GHG fluxes were measured between 21 and 23 rd July 2015, approximately 5 months after the end of felling. Three replicated soil samples were taken from 0-20 cm depth below the litter layer around each chamber at each site for bulk density, pH, total C and total N content. Soil pH was measured by mixing 5 g soil samples with 25 ml H2O with analysis using a pH meter probe . As measurements were taken from the same eight chambers through time, the analysis was conducted as a repeated measures design, with the significance of felling/non-felling differences determined against the individual chamber data (n=16), not against individual observations. A linear mixed-effects model (nlme, Pinheiro et al, 2018) was used to structure and analyse the data for all flux data. Where necessary, flux data were transformed to obtain a normal distribution. Management type (i.e. 205 felled/unfelled), temperature at 2 and 10 cm and soil moisture (plus all two-way interactions) and date (plus interaction with management type) were treated as fixed effects. Chambers were treated as a random effect (for repeated measures design). For post-felling data, residual analysis indicated that the variance was larger for chambers in the felled area and by individual chamber, therefore weighted variance structures were incorporated within the mixed effects models to account for within-type and within-chamber heterogeneity. 210 A range of autoregressive moving average corARMA (corARMA autoregressive moving average, Pinheiro et al, 2018) models were applied to each model, to account for temporal autocorrelation within chambers. Analysis of the (partial) autocorrelation function indicated potential temporal autocorrelation up to one previous time point, therefore all combinations of corARMA structure (up to one previous time point) were applied and the best fit model determined using Akaike's Information Criteria (AIC, R Core Team, 2020) applied to the maximum likelihood fits across each of the gas fluxes separately. 215 The significance of the fixed effects were subsequently determined using analysis of deviance (Anova Chi square tests, Fox and Weisberg, 2011) with non-significant interactions and effects (P>0.05) removed from the final models. (Fox and Weisberg, 2011).
Annual cumulative fluxes of CH4, N2O and CO2 were estimated to assess the inter-annual variations in each area.
Within each year, median flux values were calculated across the eight replicate chambers, along with lower and upper quartiles, 220 maximum and minimum values. These values were then accumulated across the year, by taking the mean of two consecutive flux values and multiplying it by the number of days between the measurements and summing over the monitoring period for each 12 months (as indicated in Table S1) and adjusting to 365 days. For year one where the felling interrupted the measurement in area B, the same cumulative time period was taken for both areas. Annual cumulative fluxes were converted into CO2 equivalent (CO2e) using the global warning potential (GWP) for a 100-year time horizon of 34 for CH4 and 298 for 225 N2O (IPCC, 2013).

Soil temperature, moisture
Over the 4 years, soil temperature measured during gas sampling varied between 0.1 and 26.0 o C (mean 9.7 o C) at 2 cm soil depth ( Fig. 23a) and 1.3 o and 15.5 o C (mean 8.3 o C) at 10 cm depth (data not shown), although soil temperatures were never above 15 o C under the trees in area A. Before felling in year 1, there were only small differences (p < 0.001) in mean soil 235 temperature between areas A and B (7.1 o and 8.5 o C, respectively) at 2 cm soil depth and at 10 cm (7.1 o and 7.9 o C) (Table   12), probably caused by sampling the area B later in the day than A. After felling, the mean soil temperature increased in the Commented [SY18]: Response to RC1 Table 1: we moved the  table to Supplementary information (now Table S1) felled area at 2 cm soil depth (mean 14.3 o compared to 9.2 o C in the unfelled area, p<0.001), and at 10 cm depth (10.2 o C compared to 8.4 o C in unfelled), due to removal of the tree cover (Xenakis et al., 2021). In the following two years, the soil temperature remained higher in the felled area by up to 3.4 o C compared to area A, at the 2 cm depth (1.6 o C at 10 cm) and 240 was on average 2 o C higher than before felling (there was no change in the mean soil temperature in A).  No significant differences in soil moisture content (by volume) were observed between the two areas before felling with a mean of 0.36 in area A and 0.39 m 3 m -3 in B in year 1 ( Fig. 23b and Table 12). However, after felling the soil moisture content was significantly higher in the felled area than the unfelled (mean 0.62 m 3 m -3 compared with (cf.) 0.38 m 3 m -3 , p<0.001), due to the reduced evapotranspiration after tree removal (Xenakis et al., 2021). In both areas, there was a pronounced 260 seasonal variation in soil moisture in the year before felling (2014) and in the first few months of 2015, with higher moisture in winter months than in summer, but this pattern was not clear thereafter because the rainfall was more evenly distributed during the last two years of this study (Fig. 34).

Soil parameters and root biomass measurements
Soil parameters for the 0-15 cm layer measured 5 months after felling in year 2 showed no significant differences between the unfelled and felled areas in mean soil pH (3.6 and 3.8, respectively) and soil total N content (1.8 and 1.95 %) but the mean soil total C content was about 17% lower (p <0.001) in the felled area (52.0 and 43.3 %). The soil bulk density was significantly 275 (p <0.001) higher at 0.30 g cm -3 in the felled area compared with 0.22 g cm -3 in the unfelled, but both values are typical for peaty gley soils (Vanguelova et al., 2013). The higher bulk density and lower C content could be an effect of felling, although as no pre-felling measurements were taken it may be because of site differences. Mean fine (< 2mm) live root mass was much lower in the felled area, as expected (1.6 t ha -1 cf. 4.9 t ha -1 ), but the difference was smaller for the mean fine dead root mass (0.35 and 0.53 t ha -1 , respectively), probably due to partial or complete decomposition during the 5-months after felling prior 280 to the measurements.

GHG fluxes
Large variations were observed in the fluxes between the 8 replicate chambers after felling with some high outliers, particularly for CH4 (Fig. 2c-e). Therefore, to reduce bias in the annual budget estimates of CO2, CH4 and N2O fluxes and enable a robust comparison between annual fluxes before and after felling, the annual and cumulative fluxes were based on the median of the 285 replicate chambers as described in the statistical analysis section.

CO2 effluxes
In the first year before felling, there were no significant differences in soil CO2 effluxes between areas A and B (median 1.54 and 1.75 g CO2-C m -2 d -1 , respectively, Fig. 23c). In the following 3 years after felling, CO2 effluxes became significantly (p<0.001, Table 23) lower in the felled area (median 1.10, 0.90, and 0.92 g CO2-C m -2 d -1 in year 2, 3 and 4 respectively) than 290 in the unfelled area (2.44, 1.69 and 1.44 g CO2-C m -2 d -1 ). There was a clear seasonal variation in the CO2 effluxes at both areas, which as expected followed that of soil temperature with maximum effluxes during June to September.  There was no significant correlation between CO2 effluxes and soil moisture, in part because in the felled area the moisture was high most of the time (Fig. 23b), although low CO2 effluxes were generally associated with high soil moisture 305 (and lower temperatures) in the winter and the highest effluxes were observed during low soil moisture periods in the warmer temperatures of the summer, particularly in the unfelled area. CO2 effluxes increased (p ≤ 0.002) with soil temperature at the 2 and 10 cm depths at both areas which was best described by exponential correlation relationships as shown in Fig. 45 for the periods before and after felling. Before felling, the CO2 efflux response to soil temperature was not significantly different  Figure 5. Exponential relationship between soil CO2 effluxes and soil temperature at 2 cm depth measured from area A (mature spruce stand, black) and area B (clearfell site, red), before (dashed lines, 'yr1') and after (solid lines, 'yr2-4') felling. Equations and R 2 for fitted lines shown. 320 Table 4 Table 3. Annual apparent Q10 values for soil CO2 effluxes in area A (mature spruce stand, remaining) and area B (mature spruce stand, clearfelled area after year 1) in Harwood Forest during the study period (as defined in Table S1). The annual cumulative CO2 effluxes from the felled and unfelled areas were not significantly different before felling 325 with large overlap between the 95% confidence intervals (Fig. 5a6a), and the annual CO2 effluxes were 19.8 and 24.0 t CO2e ha -1 yr −1 in areas A and B, respectively (Table 45). After felling, however, there was a clear divergence in the effluxes between the areas, with the CO2 efflux dropping sharply in the felled area (B). In the first year after felling, the annual CO2 efflux reduced to its minimum value of 8.9 compared with 23.0 t CO2e ha -1 yr −1 from the unfelled area A. In years 3 and 4 the annual effluxes in the felled area increased gradually to 11.3 and 12.2 t CO2e ha -1 yr −1 , respectively, but was still lower than the 330 unfelled area (20.3 and 18.4 t CO2e ha -1 yr −1 ).

CH4 fluxes
Soil fluxes of CH4 in both areas of the forest were generally low throughout the study period, particularly before felling with no significant differences (Fig. 3d2d). Fluxes were predominantly negative (i.e. removal from the atmosphere) in unfelled area A and before felling in area B with a median flux of -0.33 mg CH4-C m -2 d -1 and -0.21 mg CH4-C m -2 d -1 , respectively. After 345 felling, area B became a significant (p <0.001) source of CH4 and fluxes increased rapidly in the following 2 years to its maximum in year 4 (2.48 mg CH4-C m -2 d -1 ) compared to the unfelled area, which remained unchanged with a small CH4 sink (-0.33 mg CH4-C m -2 d -1 ).
Fluxes of CH4 varied more between the flux chambers after felling particularly in year 3 and 4 (Fig. 56b). Although both soil moisture and CH4 fluxes increased in the felled area after felling, the increased fluxes and the variation between 350 chambers cannot directly be related to the soil moisture as no significant overall correlation was observed. In addition, including interactions between soil moisture and temperature in the statistical model did not better explain the variation in CH4 fluxes and difference between areas (Table 32). However, the analysis showed that CH4 fluxes after felling were best explained by the soil temperature (negative association) at both the 2 cm and 10 cm depths.
Annual fluxes (expressed as GWPCO2 equivalents using the Global Warming Potential) of CH4 in the first period 355 before felling were -0.038 and -0.026 t CO2e ha -1 yr −1 from areas A and B, respectively (Table 45). In the following years, the felled area (B) became a consistent source of CH4 of 0.018, 0.194, and 0.335 t CO2e ha -1 yr −1 in year 2, 3 and 4 respectively.
In contrast, the unfelled area (A) remained a small CH4 sink with a mean GWP flux value of -0.050 t CO2e ha -1 yr −1 .

N2O fluxes
There were no significant differences in N2O fluxes between the two areas before felling (Fig. 3e2e), although the median flux 360 was higher in B (0.33 mg N2O-N m -2 d -1 and 0.55 mg N2O-N m -2 d -1 for A and B respectively). Maximum N2O fluxes of 1.30 and 2.12 mg N2O-N m -2 d -1 were measured from area A and B, respectively, between May and June 2014. After felling, N2O fluxes in the felled area B became significantly (p = 0.01) higher in year 2, 3, and 4 (median 1.83, 1.45, and 0.72 mg N2O-N m -2 d -1 ) than in the unfelled area A (0.63, 0.39, and 0.28 mg N2O-N m -2 d -1 ) with a maximum flux of 6.01 mg N2O-N m -2 d -1 measured in the first year after felling in August 2015. 365 There were no significant correlations between N2O fluxes and soil temperature before felling. However, after felling, N2O fluxes showed a seasonal pattern that followed (p <0.001) that of soil temperature at both depths in both areas with maximum fluxes during periods from June to October. Soil moisture and its interactions with soil temperature (at 10 cm) were the main driver for N2O fluxes before felling (p = 0.002 and p = 0.007, respectively). After felling, however, the soil moisture remained high (Fig. 23b) and no direct effect on N2O fluxes was observed. This could be due to the significant (p < 0.01) 370 negative correlation between soil temperature (at both depths) and moisture before felling, compared with after felling where the seasonal pattern in soil moisture was less clear (Fig. 23b and Fig. 34) and no significant correlation occurred, so that the soil temperature became the main driver of N2O fluxes.
There was a large variation between chambers in N2O fluxes throughout the study period in both areas of the forest (evident in the confidence intervals shown in Fig. 56c). Before felling, N2O annual fluxes (expressed as CO2 equivalents GWP) 375 were 0.62 and 1.03 t CO2e ha -1 yr −1 in areas A and B, respectively (Table 45), a much smaller than the GWP contribution to the total GWP from than the CO2 effluxes. After felling, the annual fluxes of N2O increased and the highest annual fluxes were measured from the felled area in the two consecutive years after felling with 2.25 and 2.52 t CO2e ha -1 yr −1 in year 2 and 3 respectively. However, at the end of the monitoring period in year 4, the annual flux of N2O returned to a rate similar to that before felling, (1.24 t CO2e ha -1 yr −1 ). 380

Effect of felling on CO2 effluxes
Felling and removal of the trees reduced soil CO2 effluxes by 55%, comparing the mean over the 3 years post-felling with that pre-felling, or 47%, comparing the clearfelled with the mature stand (Table 45). Presumably this was a consequence of the reduction in autotrophic root and rhizosphere respiration (e.g. Boone et al., 1998, Takakai et al., 2008, and measurements 385 about 5 months after felling showed a reduction from 4.9 t ha -1 to 1.6 t ha -1 in the live fine root mass (i.e. diameter <2mm).
Living fine roots and their associated mycorrhizae can contribute up to 59% of total respiration (Ewel et al., 1987;Irvine and Law, 2002;Subke et al., 2006), and for similar spruce stands in Harwood Forest about 40% contribution has been estimated previously (Zerva and Mencuccini, 2005). However, CO2 effluxes might also increase directly after felling due to an increase in decomposition of fine roots and associated ectomycorrhizal biomass and litter. Therefore, the net CO2 efflux will be 390 determined by the balance between the reduction in autotrophic root and rhizosphere respiration and the increase d decomposition (Köster et al., 20143) which can be short-lived depending on environmental factors such as soil temperature and moisture (Davidson et al., 1998;Skopp et al., 1990).
Before felling and in the unfelled stand, CO2 effluxes showed a strong seasonal pattern that followed that of soil temperature with higher effluxes during summer, as expected, when fine root density and plant growth activity are highest. 395 The apparent Q10 maximum values of these CO2 effluxes (Table 43) were lower than the value reported from a larch forest in eastern Siberia (5.92; Takakai et al., 2008), but at the higher end of those reported from a temperate deciduous oak forests in UK (1.60 to 3.92; Yamulki and Morison, 2017 and 2.2 ;Fenn et al, 2010). The response of CO2 efflux to temperature became weaker after felling (Fig. 45), even though the soil temperature was substantially higher (a decrease in the apparent Q10 values by up to 36% over the 3-year period; Table 43). This agrees with the study of Zerva and Mencuccini, (2005) who also observed 400 a weaker association of CO2 effluxes with soil temperature over a 10 month period after felling at another site in this forest.
They noted that this was probably because of the increased water content (also evident in our study, Fig. 32b), the death of fine roots and the disturbance of the soil caused by tree harvesting, and suggested that autotrophic root respiration was more responsive to temperature than heterotrophic microbial respiration. However, this effect could also be because the apparent Q10 of soil CO2 efflux determined from field measurements like these with trees present is influenced by the seasonal changes 405 in radiation and photosynthesis and their positive association with seasonal temperature, rather than an altered temperature sensitivity after felling. A larger reduction of 64% in the Q10 value of soil CO2 efflux, compared to that found in this study, was reported from a larch forest in eastern Siberia comparing a disturbed clear-cut site with a forest site (2.1 cf. 5.9, Takakai et al., 2008). At the end of our study period in year 4, the Q10 value for soil CO2 efflux in the felled area (3.46) was close to that of before felling (3.16 in B-yr1), which could be due to ground vegetation growth near but outside the gas flux chambers. 410 It may also be a result of a drop in the rainfall and soil moisture in the period between May and June in year 4 (Fig. 32b, Fig.   34) so the response to soil temperature became stronger, as in the period before felling. There may also have been a recovery from any compaction effects during harvesting, as Epron et al. (2016) showed that compaction by timber forwarding machinery after harvesting a French oak forest on a mineral soil decreased the Q10 values of soil CO2 by 16-22%. In contrast, Kulmala et al. (2014) found that the temperature dependency of the CO2 efflux was not affected by clear-cutting of a Norway spruce 415 forest on organo-mineral soil in southern Finland.
The response of CO2 efflux to soil temperature is likely to have been affected by the substantial increase in soil moisture after felling (Fig. 23b) as noted by Zerva and Mencuccini (2005) and Kulmala et al. (2014). This increase in soil moisture (with values frequently > 0.6 cm 3 cm -3 ) could have affected soil respiration by limiting the diffusion of substrates and O2 to microorganisms (Skopp et al., 1990). Yamulki and Morison (2017) could not detect an effect of soil moisture alone on 420 soil respiration in an oak forest in SE England, but the combined model of soil temperature and moisture explained 73% of the CO2 efflux variations. It is also possible that some of the CO2 produced may have dissolved in the soil water and gone undetected (Zerva and Mencuccini, 2005), but probably this effect was negligible here because of the low solubility of CO2 at the low soil pH (3.8).

Effect of felling on CH4 fluxes 425
Clearfelling changed the soil from a small annual sink of CH4 to a small net source over the 3-year monitoring period after felling, while the unfelled stand remained a sink in all years. The shift in CH4 fluxes from net uptake to net emissions by clearfelling has been reported in several other studies: Zerva and Mencuccini (2005) (2019) from Scots pine nutrient-rich peatland forest in southern Finland. As soil CH4 production requires anaerobic conditions (Conrad, 2007), this shift from sink to source was probably caused by the substantial increase in soil moisture (62% higher) and soil temperature (6 o C) in the first year after felling. An increase in soil moisture can increase the anaerobic conditions that favour CH4 production by methanogenic archea (e.g. Sundqvist et al., 2014) and therefore can change the direction of the CH4 flux. Generally, soil temperature and particularly moisture are considered to be good predictors 435 for CH4 behaviour (Lavoie et al., 2013) but disentangling the two influences is difficult in field conditions.
In this study, CH4 fluxes and the flux variations between chambers increased significantly in the following years after felling. The increase in fluxes was modest in the first year after felling but more substantial in the second and third years (Fig.   56b) which may reflect a time lag in the microbial community changing and the fungal decomposition (Glassman et al., 2018).
Although the increase in fluxes might be attributed to the substantial increase in the soil moisture, there was no direct 440 correlation between CH4 fluxes and soil moisture. Some previous studies have also attributed the increase in CH4 fluxes after felling to an increase in soil moisture or rise in the water table (Sundqvist et al., 2014;Zerva and Mencuccini, 2005;Korkiakoski et al., 2019;Epron et al., 2016) while some observed no direct correlations between CH4 fluxes and soil moisture after felling (Lavoie et al., 2013;Takakai et al., 2008;Wu et al., 2011;Mäkiranta et al., 2012;Sundqvist et al., 2014 andZerva andMencuccini, 2005). The lack of direct correlation between CH4 fluxes and soil moisture has previously been attributed to: 445 i) insensitivity of methanotrophic activity to small variations in soil moisture and temperature (Sjögersten and Wookey, 2002;Peichl et al., 2010;Wu et al., 2011;Mäkiranta et al., 2012) particularly when fluxes are relatively low such as in this study; ii) CH4 being produced at a depth greater than that of the measured soil moisture (Zerva and Mencuccini, 2005), and iii) to other overriding biological and physical factors (Lavoie et al., 2013).
Soil temperature determines CH4 fluxes by influencing methanogenic and methanotrophic activity differently (Luo 450 et al., 2013;Aronson et al., 2013). Generally, CH4 consumption by methanotrophs is less responsive to temperature than CH4 production by methanogens as consumption is mainly limited by atmospheric CH4 diffusion (Dunfield et al., 1993;Kruse et al., 1996). This is in line with the results here, as the statistical analysis showed soil temperature better explained CH4 fluxes after felling when it became a source than before when the soil was a CH4 sink. However, soil temperature can be positively Sjögersten and Wookey, 2009;Lavoie et al., 2013) depending on soil moisture and other factors that affect microbial CH4 production and consumption. We found no significant interactive effect between temperature and moisture on CH4 fluxes here.
Clearfelling can also affect other factors that play a key role in microbial CH4 production or oxidation (and CO2 and N2O fluxes) by increasing the substrate availability and soil N (NH4+ and NO3 -) concentrations (e.g. Bradford et al., 2000;Wang 460 and Ineson, 2003), for example, as a result of nitrogen release from litter or brash; or reduced soil pH (Dalal and Allen, 2008).
The higher Felling increased the bulk density in the felled area from 0.22 to (0.30 g m -3 cf. 0.22 g m -3 ) which may have been caused by compaction as a result of the machinery traffic which canmay have contributed to reduced CH4 uptake, increased CH4 production and/or CH4 release (Teepe et al., 2004;Frey et al., 2011). However, the change in soil substrates (organic matter and microorganisms) after felling was not measured in this study, and the differences between the unfelled and felled 465 area for total soil N content (1.95% and 1.83%, respectively) and soil pH (3.6 and 3.8 respectively) were small.

Effect of felling on N2O fluxes
N2O fluxes (Fig. 23e) increased after felling and there was an association with soil temperature at 10 cm depth (Table 23). Soil moisture was a significant driver for N2O fluxes before felling, but it became less significant after felling. The effect of clearfelling in increasing N2O fluxes by similar magnitudes has also been shown from forests on both mineral and peat soils 470 (e.g. Saari et al., 2009;Zerva and Mencuccini, 2005;Mäkiranta et al., 2012;Pearson et al., 2012), but an order of magnitude higher N2O flux was measured after clearfelling a nutrient-rich drained peatland forest in southern Finland (Korkiakoski et al., 2019).
N2O fluxes in temperate forest soil are generally expected to be low because of the high C:N ratios in the litter and topsoil (Butterbach-Bahl and Kiese, 2005;Jarvis et al., 2009) but can have a high spatial variability because of the variability 475 in the controlling environmental factors (Peichl et al., 2010;Fest et al., 2009). In this study, both areas of the forest showed large between-chamber variation in N2O fluxes throughout the study period. This could be due to variations in soil moisture within different flux chambers, particularly after rainfall, and the variability between chambers in soil characteristics, litter amount, mineral N availability and microbial biomass, which were not measured. N2O fluxes (Fig. 3e) increased after felling and there was an association with soil temperature at 10 cm depth (Table 3). Soil moisture was a significant driver for N2O 480 fluxes before felling, but it became less significant after felling. As the soil moisture was consistently higher after felling (Table   12), the lack of relationship with soil N2O fluxes may indicate that the moisture content was not limiting for soil N2O production by the main microbial nitrification and denitrification processes and that other factors were responsible for N2O flux variations.
As mentioned previously, the total soil N content measured some 5 months after felling was very similar to that in the unfelled area. However, microbial N2O production in the following years after felling was probably influenced by a 485 declining slow N release into the soil from the decomposition of fresh tree harvest residues and roots (Zerva and Mencuccini, 2005;Yashiro et al., 2008;Saari et al., 2009) but stimulated by warmer soil temperatures in the summer (Kulmala et al., 2014). This is consistent with the significant relationship between N2O fluxes and soil temperature, particularly at the 10 cm depth, which might imply that N2O production was more prominent at the deeper, more anaerobic soil depth, probably caused by denitrification brought about by an increased respiratory sink for O2. In the absence of plants, there is no competition for this 490 newly available N, thereby maximizing substrate availability for microbial N2O production and release (Skiba et al., 2012).
Three-fold higher N2O emissions from Finnish drained peatland pine forest plots with logging residues than from unlogged have been reported (Mäkiranta et al., 2012). Such effects of increasing soil temperature combined with microbial activities and microbial biomass N in increasing N2O fluxes have also been reported by other studies (Ineson, et al., 1991;Smolander et al.,1998;Zerva and Mencuccini, 2005;Papen and Butterbach-Bahl, 1999;Smith et al., 2018). 495 The effect of clearfelling in increasing N2O fluxes by similar magnitudes has also been shown from forests on both mineral and peat soils (e.g. Saari et al., 2009;Zerva and Mencuccini, 2005;Mäkiranta et al., 2012;Pearson et al., 2012), but an order of magnitude higher N2O flux was measured after clearfelling a nutrient-rich drained peatland forest in southern Finland (Korkiakoski et al., 2019). It is pertinent to mention here the study of Liimatainen et al. (2018) from a range of afforested northern peat soils in Finland, Sweden and Iceland, where they suggested that high N2O fluxes were linked to 500 availability of peat phosphorus (P) and copper (Cu) which could be released with other nutrients by harvesting from soil disturbance and brash (Rodgers et al., 2010), and that low P and Cu concentrations can limit N2O production even with sufficient N availability. effluxes reduced directly after felling, increasing gradually thereafter; CH4 effluxes increased sharply in year 3 and 4 after felling; and N2O effluxes increased in year 2 and 3 after felling but decreased thereafter. Currently, emissions resulting from forest management, such as felling and thinning are not considered by IPCC Tier 1 and 2 LULUCF inventories, although these 510 management activities could potentially change emission rates. In this study, the effect of clearfelling on the overall soil GHG budget and theIn order to compare the contributions of each gas to the total GHG soil flux and the effect on emissions and therefore radiative balance, we expressed the fluxes in CO2 equivalents as normally done (were assessed by comparing the annual fluxes using their global warming potential usual (GWPglobal warming potential values, ( Table 45). Before felling, the total GWP flux (sum of CO2, CH4 and N2O emissions in t CO2e ha -1 ) were was 20.4 in A and 25.0 in B areas and CO2 515 effluxes dominated the total GHG GWP, contributing up to 97%. The total GWP flux dropped in the first and second year after felling to approximately half (44% and 56% lower annual effluxes, respectively). The contribution of CO2 to the total GHG flux expressed as CO2 equivalents in the unfelled area remained constant at about 98% throughout the study period but decreased in the felled area to ca. 80% in the first and second year after felling. This was due to the doubling of N2O flux which contributed up to 20% of the total GWP flux in the two years after felling. Although the felled site became a source of 520 CH4, its contribution to the total GWP flux was always small (<2%). For the same periods, the contributions of N2O and CH4 to the total GWP flux in the unfelled area were small and remained constant at 2.9% and -0.2% respectively, similar to their values in year 1. In the last year, N2O annual efflux in the felled area halved to 1.2 t CO2e ha -1 , still 20% higher than before felling, but efflux of CH4 continued to increase to 0.34 t CO2e ha -1 . Over the 3 years since felling, the total soil GHG flux GWP was reduced by 45% (from 25.0 to 13.8 t CO2e ha -1 ) due to the much larger reduction in soil CO2 efflux than the increases in 525 N2O and CH4 fluxes. In the unfelled area there was a reduction of approximately 8%, presumably due to changing weather conditions. Fig. 7. summarises the overall results by showing the ratio of the annual soil GHG fluxes in B to that in A, with the assumption that the year-to-year variation in the unfelled area A is an indicator equally of what conditions would have been for B in the absence of felling. The figure shows a clear annual GHG flux response to felling where CO2 effluxes reduced directly after felling, increasing gradually thereafter; CH4 effluxes increased sharply in year 3 and 4 after felling; and N2O 530 effluxes increased in year 2 and 3 after felling but decreased thereafter. This study is one of very few longer-term assessments of the impacts of clearfell harvesting on the GHG balance of Sitka spruce forest. However, because of the limitations of the periodic manual closed chamber measurement technique used 540 here, it does not take into account the daily temporal flux variations or cover fluxes from ground occupied by brash or stumps.

Clearfell harvesting effect on the GHG balance
Simultaneous eddy covariance (EC) measurements (Xenakis et al., 2021) over the 3-year period after clearfelling at this study area showed approximately 3 times higher ecosystem respiration (32.8 t CO2 ha -1 yr −1 ) than our estimated mean annual soil CO2 efflux (10.8 t CO2 ha -1 yr −1 ). The difference presumably indicates high CO2 effluxes coming from the brash mats and stumps, as also suggested by Zerva and Mencuccini (2005), which could not be measured by the small chambers used in this 545 study plus above and belowground respiration of the colonising vegetation which was not included in the chambers. The differences could also be due to spatial heterogeneity over the site as the soil flux chambers are only partially sampling ground that is representative of the EC footprint. Moreover, it has been indicated that peat decomposition after felling is stimulated by nutrient release from brash (Vanguelova et al., 2010) and therefore higher N2O and CH4 fluxes may also be expected as a result of increasing mineral N release and the presence of more labile organic matter, respectively, from areas with brash. Mäkiranta 550 removal post-harvesting (e.g. for biofuel as suggested by Mäkiranta et al., 2012) might be a way of limiting GHG effluxes from peat decomposition. More information on GHG fluxes from brash and stumps, and the underlying soil processes that might be influenced by felling is a priority for future research. 555

Conclusions
In this upland Sitka spruce plantation on organic-rich peaty gley soil, clearfell harvesting affected soil GHG fluxes by increasing soil temperature and moisture content and reducing fine root mass which affected the soil nutrient and organic C supply and associated microbial populations, activities, and decomposition rates. Although soil moisture increased significantly after felling, there were no direct correlations with the soil GHG fluxes probably because there was limited 560 variation in the high soil moisture after felling. By contrast, there was a good correlation between GHG fluxes and the soil temperature which exhibited much larger temporal variation. This study does not take into account fluxes from brash decomposition, because of the small flux chamber areas, and therefore our total measured soil GHG efflux after felling probably underestimate that of the site as a whole. Over the 3-year measurement period after felling, soil CO2 effluxes reduced substantially (55%) due to cessation of root respiration outweighing increased decomposition. For the same period, CH4 fluxes 565 changed from a small net sink to a net source, increasing throughout and N2O fluxes increased substantially in the 2 years after felling. Mean soil CO2, CH4 and N2O fluxes over the 3-year period after felling contributed 83, 1 and 16%, respectively, to total GHG flux on a CO2 equivalentse basis (GWP) with an overall reduction of 45%, due to much larger soil CO2 flux reduction than the combined soil CH4 and N2O flux increases.

570
Author contributions. SY, JILM, GX and MP designed the study. GX, AA and SY carried out the gas sampling and soil environmental measurements. JB carried out the GC analyses, and SY did all the data analyses. JF was responsible for statistical analysis and SY prepared the paper, with contributions from JILM and GX.

Supplementary Information
785 Figure S1.. 2. Comparison between mean CO2 effluxes measured by gas chromatograph (GC) from the 8 static chambers in 790 each forest area (A = mature stand, B = clearfell) and mean CO2 effluxes measured by infra-red gas analyser (IRGA), from 25 static chambers in the same areas. The regression intercepts were not significantly different from zero, so the regression was through the origin. The RMSE values were 0.37 and 0.41 g CO2-C m -2 d -1 for area A and area B respectively. 795 800