Drainage Ditch Cleaning Has No Impact on the Carbon and Greenhouse Gas Balances in a Recent Forest Clear-Cut in Boreal Sweden

Ditch cleaning (DC) is increasingly applied to facilitate forest regeneration following clearcutting in Fennoscandinavia. However, its impact on the ecosystem carbon and greenhouse gas (GHG) balances is poorly understood. We conducted chamber measurements to assess the initial DC effects on carbon dioxide (CO2) and methane (CH4) fluxes in a recent forest clear-cut on wet mineral soil in boreal Sweden. Measurements were conducted in two adjacent areas over two pre-treatments (2018/19) and two years (2020/21) after conducting DC in one area. We further assessed the spatial variation of fluxes at three distances (4, 20, 40 m) from ditches. We found that DC lowered the water table level by 12 ± 2 cm (mean ± standard error) and topsoil moisture by 0.12 ± 0.01 m3 m−3. DC had a limited initial effect on the net CO2 exchange and its component fluxes. CH4 emissions were low during the dry pre-treatment years but increased particularly in the control area during the wet years of 2020/21. Distance to ditch had no consistent effects on CO2 and CH4 fluxes. Model extrapolations suggest that annual carbon emissions decreased over the four years from 6.7 ± 1.4 to 1.6± 1.6 t-C ha−1 year−1, without treatment differences. Annual CH4 emissions contributed <2.5% to the carbon balance but constituted 39% of the GHG balance in the control area during 2021. Overall, our study suggests that DC modified the internal carbon cycling but without significant impact on the carbon and GHG balances.


Introduction
The boreal forest biome is an important global reservoir of terrestrial carbon (367-1716 Pg C) [1]. However, its biogeochemical cycling and thus carbon storage is sensitive to forest management practices [2]. Specifically, tree growth is often hampered by waterlogging in humid soils, which is also a typical feature in boreal upland forest areas [3][4][5]. Therefore, artificial drainage ditching has been performed during the past century in large parts of northern Europe (Finland, Sweden, Estonia and Latvia) to increase tree biomass production [6]. Most of these earlier ditched forests are now reaching the end of the rotation period, in which time, the drainage capacity and efficiency in lowering the water table level (WTL) has deteriorated for many of these old ditches. In addition, decreased evapotranspiration following clear-cutting further supports an increase in the WTL that limits the establishment, survival and growth of subsequent seedling generation [7,8].
As one solution for mitigating increasing soil wetness, ditch cleaning (DC) following forest harvest may help to restore the drainage function of the ditches and thereby regain desired tree growth rates (Paavilainen and Päivänen, 1995) [3]. According to the official statistics from the Finnish National Forest Programme and the Swedish National Forest Inventory (NFI), about 65,000 ha and 10,000 ha of forest lands have been ditch cleaned (1) Quantify the magnitudes of CO 2 and CH 4 fluxes from seasonal to inter-annual scales; (2) Investigate the effects of DC on the spatio-temporal variations in CO 2 and CH 4 fluxes; (3) Identify environmental factors that drive the changes in CO 2 and CH 4 fluxes in response to DC; (4) Estimate the effect of DC on the annual C and GHG balances.

Site Description and Experimental Design
This study was conducted at a forest clear-cut site ("the Pettersson site" 64 • 12 56 N, 20 • 49 32 E , 60 m.a.s.l.) located approximately 3 km north of Robertsfors town in the county of Västerbotten, Sweden. The local climate according to the Köppen-Geiger classification [34] is characterized as boreal (Dfc; also called subarctic) with persistent snow cover during~6 months. According to the data from the nearest weather stations (Brände and Bjuröklubb, 15 km and 46 km from the site, respectively) of the Swedish Meteorological and Hydrological Institute (SMHI; www.smhi.se, accessed on 1 February 2022), the 30-year (1991-2020) mean annual temperature (Recorded in Bjuröklubb) is 4.0 • C, and the mean annual precipitation (Recorded in Brände) is 701 mm, of which about 35% falls as snow.
While detailed Records for the timing of the original drainage and plantation of the former stand are missing, historical air photos suggest that trees and the ditch network of the study site already existed during the early 1960s. Prior to harvest in October 2016, the former stand had a tree volume of 274 m 3 ha −1 and consisted of about 75% pine (Pinus sylvestris L.), 22% Norway spruce (Picea abies L.) and 3% broad-leaved tree species (e.g., Betula spp). Soil preparation by mounding was conducted across the entire site in August 2018, followed by the establishment of tree seedlings in June 2019, which consisted of 30% Norway spruce and 70% pine, planted with density of 2100 plants per hectare.
The drainage ditch network is composed of two main ditches that diverge water into different directions, i.e., south-east and south-west ( Figure 1). In January 2020, DC was performed along the ditch network branch draining towards the south-east (its surrounding area is hereafter referred to as the DC area), using a tracked excavator. During DC, both vegetation and deposited material were removed from the ditches and piled up along both sides of the ditches. After DC, the open ditches had a trapezoidal shape, with a depth of about 1 m and a width of about 0.5 at the bottom and towards 2 m at the top. The ditch network branch draining to the south-western side was left uncleaned (hereafter referred to as the control area) and characterised by stagnant surface waters due to sediment deposition and vegetation in-growth, which reduced its drainage functions.

Greenhouse Gas Flux Measurements
We conducted CO2 and CH4 flux measurements using the closed dynamic chamber method [35]. The measurements were carried out during daytime in approximately biweekly intervals during the snow-free period (May to October) from 2018 to 2021. In May 2018, square aluminium frames (48.5 × 48.5 cm) with a frame base down to 5 cm below the soil surface were permanently installed at each plot. The observed vegetation species and coverage within the frames were consistent with the surrounding area indicated in Section 2.1 over the study years.
Forest-floor CO2 and CH4 fluxes were determined by placing a chamber connected in a closed loop to a portable GHG analyser (Gas Scouter G4301, Picarro Inc., Santa Clara, CA, USA) onto the pre-installed frames. The Gas Scouter has a built-in sampling pump circulating air between chamber and analyser at a flow rate of 0.001 m 3 min −1 , and Records CO2 and CH4 concentrations at a frequency of approximately ~1 Hz, with precision levels of ±150 ppb and ±0.8 ppb, respectively. In the first two years (2018/19), we used a transparent chamber with the dimensions of 48.5 × 48.5 × 30 cm. During the subsequent two years (2020/21), a taller transparent chamber (48.5 × 48.5 × 50 cm) was required to cover the growing herbaceous vegetation within the frames. While effective mixing of the cham- The site is located in a region where Podzol is the dominant soil type (Geological Survey of Sweden (SGU)), with soil layers consisting of thin incoherent surface layers of postglacial sand and gravel, while finer soil particles of clay and silt are dominant at 0.5 m below the ground surface. The thin organic layer at our study site was partly mixed into the upper mineral soil layer by the harvest machinery, which resulted in elevated carbon and nitrogen concentrations of 25.1 ± 5.3 % and 0.84 ± 0.19 % (i.e., C:N ratio: 31 ± 1.5) in the 0 to 10 cm profile, respectively. In comparison, carbon and nitrogen content decreased to 10.6 ± 4.5 % and 0.40 ± 0.17 % (i.e., C:N ratio: 31 ± 1.1), respectively, within 10 and 20 cm depth. There was no significant difference in carbon and nitrogen concentrations between the control and DC areas based on the soil samples taken from 24 measurement plots in August 2018. The main herbaceous plant species that were established across the site during the study years following harvest included Deschampsia flexuosa L., Trichopho-Forests 2022, 13, 842 4 of 24 rum alpinum L., Epilobium angustifolium L., Potentilla erecta L., Vaccinium myrtillus L. and Vaccinium vitis-idaea L.
The experimental design included two parallel transects (about 30 m apart) in each of the two treatment areas with sampling plots at three distances (4, 20, and 40 m) on both sides of the ditch along each transect (i.e., 2 treatment areas × 2 transects × 3 distances × 2 sides = 24 plots) ( Figure 1). The 40 m plots were also at similar distance to the surrounding intersecting ditches (Figure 1). Pre-treatment data were collected for two years (2018-2019) followed by two years (2020 and 2021) of measurements after the ditch was cleaned in the designated DC treatment area. This experimental setup allowed us to control for potential confounding effects from both spatial (e.g., topography) and temporal (e.g., weather patterns) factors.

Greenhouse Gas Flux Measurements
We conducted CO 2 and CH 4 flux measurements using the closed dynamic chamber method [35]. The measurements were carried out during daytime in approximately biweekly intervals during the snow-free period (May to October) from 2018 to 2021. In May 2018, square aluminium frames (48.5 × 48.5 cm) with a frame base down to 5 cm below the soil surface were permanently installed at each plot. The observed vegetation species and coverage within the frames were consistent with the surrounding area indicated in Section 2.1 over the study years.
Forest-floor CO 2 and CH 4 fluxes were determined by placing a chamber connected in a closed loop to a portable GHG analyser (Gas Scouter G4301, Picarro Inc., Santa Clara, CA, USA) onto the pre-installed frames. The Gas Scouter has a built-in sampling pump circulating air between chamber and analyser at a flow rate of 0.001 m 3 min −1 , and Records CO 2 and CH 4 concentrations at a frequency of approximately~1 Hz, with precision levels of ±150 ppb and ±0.8 ppb, respectively. In the first two years (2018/19), we used a transparent chamber with the dimensions of 48.5 × 48.5 × 30 cm. During the subsequent two years (2020/21), a taller transparent chamber (48.5 × 48.5 × 50 cm) was required to cover the growing herbaceous vegetation within the frames. While effective mixing of the chamber headspace was created by the continuous recirculation pump from the analyser in the smaller chamber used in 2018/19, the larger chamber was equipped with a small fan to further support mixing in the larger and more densely vegetated headspace. A comparison of repeated (i.e., within <5 min) fluxes measured with the two different chamber sizes suggested good agreement for CO 2 (r 2 = 0.95) and CH 4 (r 2 = 0.58) fluxes, respectively. The lower r 2 for the CH 4 flux was likely due to the disturbance from the initial measurement, which may have modified the comparatively small CH 4 concentration gradient before the subsequent measurements. The mixing ratios of CO 2 and CH 4 are reported as dry air mole fractions.
During each flux measurement, a transparent chamber was used at first to measure the net ecosystem exchanges for CO 2 (i.e., NEE) and CH 4 during 90-120 s under natural light conditions. After 2 min of venting, the chamber was placed onto the same frame again and covered by an opaque blanket to estimate the CO 2 flux during dark conditions, i.e., ecosystem respiration (Reco). The gross primary productivity (GPP) was then calculated using the difference between NEE and R eco.
The average rate of the change in the mixing ratios over time (dC/dt; ppm s −1 ) was estimated using simple linear regression over a chosen data range. The flux rate was then calculated based on dC/dt as a function of chamber headspace volume, air temperature and pressure according to the ideal gas law. Poor-quality flux data were eliminated under the criteria of the root-mean-square error (RMSE) and r 2 of the chosen dC/dt. Based on visual examination of the data, CO 2 fluxes with RMSE > 2.5 ppm and r 2 < 0.90, and CH 4 fluxes with RMSE > 2.5 ppb and r 2 < 0.90 were removed. These quality control procedures led to the removal of about 2% and 5% of all CO 2 and CH 4 fluxes measured, respectively. The sign convention in this study is such that positive and negative flux values refer to a net source and sink of the GHG, respectively. We also conducted measurements of N 2 O fluxes during early and late August 2020 with static chambers following the methodology described in Appendix A. The results suggested negligible N 2 O fluxes of 14 ± 8 µg-N m −2 h −1 and 18 ± 5 µg-N m −2 h −1 in the control and DC area, respectively, which accounted for <1.1% of the GHG balance even when accounting for the 298 higher warming potential of N 2 O relative to CO 2 . We therefore did not include these measurements into the ordinary measurement program.

Measurements of Abiotic Factors
To investigate the impact of DC on the environmental factors and to explore their relationships with the measured GHG fluxes, we Recorded a suite of environmental variables in parallel during each flux sampling campaign. Specifically, manual WTL measurements were made inside PVC groundwater tubes (Ø = 32 mm external and 26 mm inside, h = 125 cm, 3 mm holes every 2.5 cm) inserted to a depth of 1 m adjacent to each flux measurement frame. Soil moisture within the upper 5 cm was measured at three sides around the frame using a GS3 combined moisture-temperature sensor connected to a handheld data logger (ProCheck, Decagon Devices, Pullman, WA, USA). Soil temperature (T s ) was also Recorded manually next to each frame at 5 and 10 cm depths (T s5 , T s10 ). Air temperature (T a ) was measured using two Hobo ® pendant temperature loggers (Onset Computers, Bourne, MA, USA). One was placed at 2 m height above surface and shaded from direct sunlight, Recording data at 15 min intervals continuously over the year, whereas the other one was operated in the chamber headspace at 5 s intervals during chamber closures. Photosynthetically active radiation (PAR) was continuously measured using two Hobo ® pendant radiation sensors with built-in loggers (Onset Computers, Bourne, MA, USA). Specifically, one sensor Recorded PAR at 5 s intervals inside the chamber during the period of a flux measurement, while another sensor was situated at the centre of the study site to log ambient PAR all year round at an hourly interval. Soil temperature at 5 cm and 10 cm depths was also continuously measured with automated TOMST ® TMS-4 sensors (TOMST, Prague, Czech Republic) installed at 6 selected plots, with each located 4 m, 20 m and 40 m from the ditches in both the control and DC areas. The hourly ambient T a and PAR data served as input for the nonlinear regression models to predict hourly CO 2 fluxes (see Section 2.4).

Vegetation Characteristics
The growth of the field-layer vegetation within each flux measurement frame was monitored during the measurement years by taking overhead images during each flux measurement campaign. These images were then analysed to derive a vegetation greenness index defined by the green chromatic coordinate (g cc ) [36][37][38] (Equation (1)).
g cc refers to the greenness index from the image taken on the frame; R, G and B denote the intensity (0-255) of the red, green and blue image channels. The g cc values were averaged for each image pixel located within the chamber frame. Based on the overhead images, we also calculated the ground vegetation areal coverage, which is defined as the vertical projection of vegetation onto a unit of land.
In July 2019 and July 2021, the ground vegetation g cc and areal coverage for the area surrounding the frames were studied by taking overhead images at 12 spots evenly distributed within 15 m radii around each flux measurement frame. The purposes of these measurements were to (1) evaluate if the vegetation cover inside the flux measurement frames was representative of that within the surrounding area and (2) to increase the sample size to test for DC impacts on ground vegetation development.
Continuous measurements of vegetation growth were also made using two phenocams (TimeLapseCam, Wingscapes, Calera, AL, USA) which Recorded images at hourly intervals with a nadir angle of 15 • towards a northernly direction. The g cc and areal coverage calculated from these images were used to calibrate and interpolate the manual overhead images to frame-specific continuous phenology time series. The latter were then used as model input for estimating annual CO 2 flux budgets (see Section 2.4).

Statistical Analysis
We first applied mixed effect models with repeated measures to test for the significant effects of the DC and distance to ditch treatments on the spatial variation of environmental (i.e., T a, T s10 , T s5 , SM, WTL, g cc and vegetation areal coverage) and flux (i.e., GPP, Reco, NEE and CH 4 ) variables ( Figure 2). The statistical models applied were as follows (Equation (2)): where y ij denotes the environmental or gas flux variable for sampling occasion i with DC treatment j (j = control or DC); β 0 denotes the overall mean of the environmental or gas flux variable; T j denotes the fixed effect of DC treatment j; β 1 denotes the fixed effect of the sensitivity to distance to ditch d; S ij denotes the random effect of sampling occasion i; ε ij denotes the random error for sampling occasion i with treatment j. The model takes into account the random effects presenting a covariance structure where correlations decrease with time [39]. Mixed effect models were proven to be robust to various data distributions [40]. We set the statistical error levels to α = 0.05 for the mixed effect models, and the standard error (±SE) of the sample means was used to indicate the level of uncertainty.

Modelling of Annual CO2 and CH4 Flux Budgets
Nonlinear regression models following [36,[44][45][46] were used to predict hourly and annual Reco and GPP fluxes based on data of Ta, PAR and gcc. Particularly, GPP from each frame was fitted to hourly PAR and frame-specific gcc data using a hyperbolic PAR function adjusted by normalised frame-specific representing seasonal variations in vegetation biomass (Equation (3)): where GPP denotes the hourly gross primary production (mg m −2 h −1 of CO2-C); α denotes the initial slope of quantum use efficiency of photosynthesis (mg µmol photons −1 of CO2-C); PAR denotes the mean of hourly photosynthetically active radiation (µmol m −2 s −1 ); Pmax denotes the maximum photosynthesis under infinite PAR (mg m −2 h −1 of CO2-C); and is the daily frame-specific chromatic greenness index (gcc) surrounding the frames normalised to the scale of 0 to 1. Next, since the mixed effect models did not provide information on the associations between environmental and flux variables, we complemented our statistical analysis by conducting a principal component analysis (PCA) with the goal to identify the three-way interactions between the treatments (DC and distance), environmental factors and flux variables (CO 2 and CH 4 ) ( Figure 2). Principal components (PCs) were calculated using the growing season means of the variables from each measurement plot. Separate PCAs were carried out for the datasets from pre-DC and post-DC treatment periods. Variables were normalised by subtracting each value from the mean and dividing by the standard deviation before inputting them to the PCA models [41]. Significant PCs were selected and presented using the Kaiser criterion [42]. The correlation coefficients (loadings) of the variables with the significant PCs were compared to identify the relationships among variables [43]. All statistical analyses were conducted using the Mathworks Matlab software R2021a.

Modelling of Annual CO 2 and CH 4 Flux Budgets
Nonlinear regression models following [36,[44][45][46] were used to predict hourly and annual R eco and GPP fluxes based on data of T a , PAR and g cc . Particularly, GPP from each frame was fitted to hourly PAR and frame-specific g cc data using a hyperbolic PAR function adjusted by normalised frame-specific g cc representing seasonal variations in vegetation biomass (Equation (3)): where GPP denotes the hourly gross primary production (mg m −2 h −1 of CO 2 -C); α denotes the initial slope of quantum use efficiency of photosynthesis (mg µmol photons −1 of CO 2 -C); PAR denotes the mean of hourly photosynthetically active radiation (µmol m −2 s −1 ); P max denotes the maximum photosynthesis under infinite PAR (mg m −2 h −1 of CO 2 -C); and g cc norm is the daily frame-specific chromatic greenness index (g cc ) surrounding the frames normalised to the scale of 0 to 1. Hourly estimates of R eco were modelled using an exponential function of soil temperature based on [47] adjusted by the normalised frame-specific g cc (Equation (4)): where R eco denotes hourly ecosystem respiration (mg m −2 h −1 of CO 2 -C); R 0 denotes soil respiration at 0 • C (mg m −2 h −1 of CO 2 -C); T s5 denotes soil temperature at 5 cm depth ( • C); b denotes the sensitivity of R eco to T s5 ; and β is a scaling parameter associated with the contribution of autotrophic respiration (R) a to R eco. Hourly estimates of GPP and R eco were then summed up for the entire year, and annual NEE was estimated by the difference of the modelled annual GPP and R eco.
To increase the robustness of the models, we pooled the available data to develop models based on the two pre-DC (2018/19) and two post-DC (2020/21) years separated for the control and DC areas, respectively. The model parameters of the final models selected based on the highest R 2 are summarised in Table 1. It is noteworthy that despite similar greenness and areal cover (Table S1), the P max in the GPP model and the β in the R eco model were higher in the DC plots during the post-DC period compared to the control plots. This is explained by the greater vascular biomass noted in the DC plots (based on visual observation) at which the same greenness and areal coverage facilitate higher photosynthetic rate compared to moss-dominated ground vegetation. During the snow period, R eco was represented by R o , which is an iterated parameter in the respiration model that describes the respiration rate when soil temperature was at 0 • C. The modelled GPP values were zero in response to the negligible exposure of ground vegetation during snow coverage (i.e., g cc norm = 0; Equation (3)).
The uncertainty of the annual CO 2 flux budgets derived by the model extrapolations was evaluated using Monte Carlo simulations. For that purpose, we assigned a normal distribution to each model input parameter (see Table 1) in accordance with its mean and standard deviation derived during model development [48]. Then, a large number (1000) of random draws were taken from the normal distributions of each model parameter. The standard deviation for the set of 1000 predicted estimates was then used to define the uncertainty of the annual CO 2 flux budgets. Table 1. Model parameters for estimating gross primary production (GPP) (Equation (3)) and ecosystem respiration (R eco ) (Equation (4)) for control and ditch cleaning (DC) clear-cut area, applied separately in the pre-DC years (2018/19) and post-DC years (2020/21); α is the initial slope of the lightuse photosynthetic efficiency (mg µmol −1 photons of CO 2 -C); P max is the maximum photosynthetic rate at light saturation (mg m −2 h −1 of CO 2 -C); R 0 is respiration rate (mg m −2 h −1 of CO 2 -C) at 0 • C; b is the sensitivity of R eco to soil temperature (T s ); β represents the contribution of autotrophic respiration to R eco ; numbers in parentheses indicate standard error; R 2 denotes the coefficient of determination of the model. Number of observations ranges between 186 and 272 for the development of each model.

Time
Pre Due to the weak relationship between CH 4 fluxes and environmental variables, annual CH 4 balances were interpolated from the median of measured CH 4 fluxes. We propagated the standard errors to estimate the uncertainty of the annual CH 4 budgets due to spatial variability among sampling plots. A baseline winter CH 4 emission rate was estimated based on the assumption that the low CH 4 fluxes Recorded from the late October measurement campaigns (see Section 3.2) continued during the following winter period (i.e., November-April).

Environmental Data
The annual mean air temperature at the nearby SMHI weather station in Bjuröklubb (46 km from the study site) was close to the 30-year average (4.0 • C in 1991-2020) in 2018 (4.4 • C), 2019 (3.9 • C) and 2021 (4.1 • C), but considerably warmer in 2020 (5.9 • C) due to an unusually warm winter period (−0.4 • C and −2.1 • C in January and February, with reference to 30-year average of −5.1 • C and −6 • C). The annual total precipitation in 2018 and 2019 at the nearby SMHI weather station in Brände (15 km from the study site) was 505 and 652 mm, which was lower than the 30-year long-term average (1991-2020) of 701 mm. In the post-DC years of 2020 and 2021, the annual total precipitation was 958 and 1007 mm, substantially higher than the long-term average. T a Recorded during the flux measurement campaigns ranged from a minimum of 5.9 ± 0.4 • C during the beginning or the end of the growing season to a maximum of 29.1 ± 1.3 • C in June or July during the four years ( Figure 3a). T s5 and T s10 measured during flux measurements ranged from 4.4 ± 0.3 • C to 19.3 ± 0.9 • C and from 4.7 ± 0.4 • C to 16.2 ± 0.7 • C, respectively, for the four years ( Figure 3a). The difference in T s5 and T s10 between the control and DC areas was <0.4 • C during the study years. T s5 and T s10 at 4 m distance from the ditch were 0.2 to 0.9 • C lower than at 40 m from the uncleaned ditches in the control area (Table S1). The magnitudes and temporal dynamics of PAR were similar during the study years ( Figure 3b). °C in June or July during the four years (Figure 3a). Ts5 and Ts10 measured during flux measurements ranged from 4.4 ± 0.3 °C to 19.3 ± 0.9 °C and from 4.7 ± 0.4°C to 16.2 ± 0.7 °C, respectively, for the four years ( Figure 3a). The difference in Ts5 and Ts10 between the control and DC areas was <0.4 °C during the study years. Ts5 and Ts10 at 4 m distance from the ditch were 0.2 to 0.9 °C lower than at 40 m from the uncleaned ditches in the control area (Table S1). The magnitudes and temporal dynamics of PAR were similar during the study years (Figure 3b).  There was no clear seasonal pattern on the WTL fluctuation, but a shallower WTL was usually observed at the beginning or the end of the growing season. In 2018, the mean WTL averaged over all flux plots and sampling campaigns was −39 ± 2 cm and −43 ± 3 cm in the control and DC areas, respectively ( Figure 3c). Thus, the mean WTL was already 5 cm higher in the control area than in the DC area during the pre-DC period. In 2019, the mean WTL rose to −27 ± 1 cm in the control area and −32 ± 1 cm in the DC area, suggesting a significant (p < 0.05) but similar increase in both areas during the second pre-DC treatment year (Table S1). After DC, the mean WTL in the control area rose to −19 ± 1 cm in 2020 and −15 ± 1 cm in 2021, whereas the WTL in the DC area remained at a similar level of −37 ± 2 cm in 2020 and −32 ± 1 cm in 2021, relative to the pre-DC year of 2019. Thus, following DC, the WTL in the DC area was 18 ± 2 cm and 17 ± 2 cm lower than in the control area in 2020 and 2021, respectively. This implies that the WTL difference between the two areas increased by 12 ± 2 cm due to the DC effect. In the control area, significantly shallower WTL was observed in the plots closer to the ditch during all four study years. In the DC area, however, plots closer to the ditch were already characterised by deeper WTL in the pre-DC years, and this spatial difference remained significant in the post-DC years (Table S1).
Similar spatial and temporal patterns were observed for the soil moisture content. Specifically, in the first two pre-treatment years, the annual mean soil moisture remained in the range of 0.29 to 0.36 m 3 m −3 in the two treatment areas (Figure 3d; Table S1). The difference between the control and DC area increased substantially in the two post-DC years, where in the control area, the soil moisture increased to 0.44 ± 0.02 and 0.43 ± 0.01 m 3 m −3 in the year 2020 and 2021. In the DC area, the soil moisture remained at 0.30 ± 0.01 and 0.26 ± 0.01 m 3 m −3 in the year 2020 and 2021, respectively, and was not statistically different from the pre-DC period. Altogether, the DC effect on soil moisture content was estimated at 0.12 m 3 m −3 based on the additional differences noted between the control and DC areas after DC. The effect of distance to ditches on the soil moisture was similar to its effect on the WTL, with soil moisture decreasing with distance to ditches in the control area but increasing with distance to ditches in the DC area during all measurement years.
In July 2019 and 2021, the seasonal maxima of the greenness index and vegetation areal coverage were not significantly different between the two treatment areas (Figure 4). Averaged over all frames and treatments, the greenness index increased from 0.24 ± 0.02 in July 2019 to 0.53 ± 0.04 in July 2021, while the areal coverage increased from 0.21 ± 0.04 in July 2019 to 0.54 ± 0.05 in July 2021. In the surrounding area, the maximum greenness index and areal coverage over all treatment plots increased from 0.22 ± 0.01 to 0.46 ± 0.02 and from 0.20 ± 0.02 to 0.57 ± 0.02, respectively, over the two years. The greenness index and vegetation areal coverage in the surrounding area at the 4 m distance from the ditch locations were, however, 15% and 30% greater in the control area than in the DC area in July of 2019 and 2021, whereas no difference between both areas was noted at the 20 m and 40 m plots (Table S1). Furthermore, the greenness index and vegetation areal coverage of frames and their surrounding areas agreed well overall, except the greenness index in the frame was 16% greater (p = 0.047) than in the surrounding areas in the DC area in July 2021. In the model extrapolation for annual gas balances, this difference was accounted for by developing treatment-specific models (see Section 2.4).

Temporal Variations of CO 2 and CH 4 Fluxes
The magnitudes and seasonal variations of daytime net CO 2 exchange (NEE) and its components (GPP and R eco ) increased throughout the four measurement years in both the control and DC areas. When averaged over all plots, the instantaneous daytime NEE switched from net emissions of 72 ± 9 mg C m −2 h −1 in 2018 to a net uptake of 168 ± 27 mg C m −2 h −1 in 2021 (Figure 5a; Table 2). From 2018 to 2021, the peak growing season GPP increased about 13 times in the DC area, relative to a 10-time increase in the control area (Table 2). Similarly, daytime measurements of R eco increased by more than three times in the DC area from 2018 to 2021, whereas in the control area, the increase was by about two times. The boxes denote the interquartile range with the median indicated inside each box, whereas the lines extending from the boxes indicate the range of each category. Different letters above boxes (a, b, c) refer to compact letter display of one-way ANOVA that indicate significant differences (α = 0.05) among groups , followed by the least significant difference test.

Temporal Variations of CO2 and CH4 Fluxes
The magnitudes and seasonal variations of daytime net CO2 exchange (NEE) and its components (GPP and Reco) increased throughout the four measurement years in both the control and DC areas. When averaged over all plots, the instantaneous daytime NEE switched from net emissions of 72 ± 9 mg C m −2 h −1 in 2018 to a net uptake of 168 ± 27 mg C m −2 h −1 in 2021 (Figure 5a; Table 2). From 2018 to 2021, the peak growing season GPP increased about 13 times in the DC area, relative to a 10-time increase in the control area (Table 2). Similarly, daytime measurements of Reco increased by more than three times in the DC area from 2018 to 2021, whereas in the control area, the increase was by about two times. The boxes denote the interquartile range with the median indicated inside each box, whereas the lines extending from the boxes indicate the range of each category. Different letters above boxes (a, b, c) refer to compact letter display of one-way ANOVA that indicate significant differences (α = 0.05) among groups, followed by the least significant difference test.
In the two pre-treatment years, CH 4 fluxes varied within the 10-90 percentile range of −0.06 to +0.07 mg C m −2 h −1 across all plots, with a median of −0.02 mg C m −2 s −1 , indicating the majority of fluxes were negative (i.e., uptake) at the site ( Figure 6). Occasional emission spikes were observed at individual plots across both the experimental areas throughout the four years ( Figure 6). In the two years after DC, the number and magnitude of CH 4 emission spikes increased, particularly during the peak season in both the control and DC areas, but with considerably higher ranges and frequencies observed in the control area. Specifically, the 10-90 percentile ranged from −0.07 to 3.6 mg C m −2 s −1 in the control area and from −0.14 to 0.11 mg C m −2 s −1 in the DC area during the two years after DC. The median CH 4 flux was 0.08 mg C m −2 s −1 in the control area and −0.05 mg C m −2 s −1 in the DC area (Table 2). It is further noteworthy that most emission spikes were Recorded at the plots 4 m from the uncleaned ditch in the control area, with median emissions of 1.2 mg C m −2 s −1 ( Table 2).

DC Effect on CO 2 and CH 4 Fluxes
During the pre-treatment period, the mixed effect models suggested that the growing season means of daytime NEE exhibited no significant differences (p > 0.05) between the control and DC areas. The growing season means of the measured GPP and R eco fluxes also did not differ significantly (p > 0.05) between the two treatment areas during the two pre-treatment years. After DC, both measured R eco and GPP became significantly (p < 0.01) greater as compared to the control area. Due to similar changes in R eco and GPP, however, NEE was not significantly affected by DC (p > 0.05).
The mixed effect models indicate no significant difference (p > 0.05) in CH 4 fluxes between the control and DC areas before the treatment ( Table 2). In the two years after DC, emissions from the control area became significantly higher (p < 0.01) than from the DC area, consistent with the increased emission spikes in the control area.  the control area. Specifically, the 10-90 percentile ranged from −0.07 to 3.6 mg C m −2 s −1 in the control area and from −0.14 to 0.11 mg C m −2 s −1 in the DC area during the two years after DC. The median CH4 flux was 0.08 mg C m −2 s −1 in the control area and −0.05 mg C m −2 s −1 in the DC area (Table 2). It is further noteworthy that most emission spikes were Recorded at the plots 4 m from the uncleaned ditch in the control area, with median emissions of 1.2 mg C m −2 s −1 ( Table 2).

DC Effect on CO2 and CH4 Fluxes
During the pre-treatment period, the mixed effect models suggested that the growing season means of daytime NEE exhibited no significant differences (p > 0.05) between the control and DC areas. The growing season means of the measured GPP and Reco fluxes also did not differ significantly (p > 0.05) between the two treatment areas during the two pretreatment years. After DC, both measured Reco and GPP became significantly (p < 0.01) greater as compared to the control area. Due to similar changes in Reco and GPP, however, NEE was not significantly affected by DC (p > 0.05).
The mixed effect models indicate no significant difference (p > 0.05) in CH4 fluxes between the control and DC areas before the treatment ( Table 2). In the two years after

Distance to Ditch Effect Effects on CO 2 and CH 4 Fluxes
The distance to ditches had contrasting effects on NEE, R eco and GPP. Specifically, NEE shifted from a net CO 2 sink to a net CO 2 source with increasing distance to ditches in the control area, whereas the magnitude of net CO 2 uptake decreased with distance to ditches in the DC area (Table 2). This spatial pattern persisted throughout all four measurement years. Meanwhile, the significant effects of distance to ditches on R eco and GPP were only noted in the control area, suggesting a decrease in their magnitude with increasing distance to the ditches. In the DC area, however, the highest R eco and GPP values were observed at 20 m from the ditch ( Table 2). These distance effects were consistent throughout all measurement years.
As the CH 4 emission spikes were concentrated at the plots 4 m from the uncleaned ditches, the mixed effect models also suggest a negative effect of distance to ditch on the CH 4 fluxes. However, additional inspection of the data suggests that this distance to ditch effect was mainly driven by a strong gradient in the control area, where CH 4 emissions, coinciding with highest WTL, remained highest at the plots closest (4 m) to ditches.

Three-Way Interaction between Ditch Treatments, Environmental Factors and GHG Fluxes
The PCA on all variables Recorded during the pre-DC years showed that in total, 83% of the total variance was explained by five significant PCs (Table S3). Together, the first two most significant PCs explained 52% of the total variance (Figure 7). PC1 revealed that vegetation growth (g cc and areal coverage) was most strongly coherent with the CO 2 component fluxes of R eco and GPP. Vegetation variables were also positively related to soil carbon and nitrogen content and negatively related to soil temperature. PCA results further showed that R eco and GPP were negatively associated with the distance to ditch and CN ratios. CH 4 flux had a strong positive connection with WTL and soil moisture in both PC1 and 2 but was not associated with the distance to ditch. During the pre-DC years, the DC treatment variable had relatively small variable loadings in both PC1 and 2, implying that pre-DC environmental conditions and fluxes were comparable between the two measurement areas.
The PCA on all variables Recorded during the pre-DC years showed that in total, 83% of the total variance was explained by five significant PCs (Table S3). Together, the first two most significant PCs explained 52% of the total variance (Figure 7). PC1 revealed that vegetation growth (gcc and areal coverage) was most strongly coherent with the CO2 component fluxes of Reco and GPP. Vegetation variables were also positively related to soil carbon and nitrogen content and negatively related to soil temperature. PCA results further showed that Reco and GPP were negatively associated with the distance to ditch and CN ratios. CH4 flux had a strong positive connection with WTL and soil moisture in both PC1 and 2 but was not associated with the distance to ditch. During the pre-DC years, the DC treatment variable had relatively small variable loadings in both PC1 and 2, implying that pre-DC environmental conditions and fluxes were comparable between the two measurement areas.  Table S3. Filled symbols denote the component loadings of the measured variables. Crosses highlighted in black (control area) and red (DC area) denote the component scores of each measurement plot. Colour intensity of component scores decreases with the distance to the ditches. Note that ecosystem uptake in NEE and GPP is given a negative sign, resulting in negative correlation to both PC1s despite positive causal relation. Abbreviations represent greenness index (g cc ), soil moisture (SM), air temperature (T a ), soil 5 or 10 cm depth temperature (T s5 , T s10 ), vegetation areal coverage (VC), water table level (WTL), soil carbon (C) and nitrogen content (N) along with their C:N ratio (CN) for environmental variables (black symbols); net ecosystem exchange (NEE), ecosystem respiration (R eco ), gross primary productivity (GPP) and methane (CH 4 ) for flux variables (blue symbols); and DC treatment and distance to ditch variable (purple symbols). Note that the DC variable was fitted as a binary variable for control (0) and DC (1) treatments.
During the post-DC years, 80% of the total variance was explained by four significant PCs (Table S3), while together, the first two most significant principal components explained 54% of the total variance (Figure 7). Compared to the pre-DC years, the DC treatment effect had a larger association with the environmental and flux variables, as reflected by the increased variable loadings in the first two PCs. Specifically, the DC treatment was associated with lower soil moisture, lower WTLs and with increased soil temperature, in correspondence with lower CH 4 fluxes. In contrast, the CO 2 component fluxes were independent from the DC effects on soil temperature and moisture levels but instead were primarily controlled by vegetation growth. The latter also remained positively related to soil carbon and nitrogen contents. The results further show that soil moisture, WTL and thereby CH 4 emissions decreased with greater distance from the ditches. Meanwhile, the small loading of distance to ditch indicated its limited association with the other environmental and CO 2 flux variables relative to the pre-DC period.

Total Annual Carbon and GHG Balances
Our model estimates suggest that our study site was an annual source of carbon during the four study years, regardless of the DC treatment ( Figure 8; Table S2). However, the source magnitude decreased by 78% and 74% in the control and DC areas, respectively, from 6.1 ± 1.2 t-C ha −1 year −1 and 7.2 ± 1.6 t-C ha −1 year −1 in 2018 to 13.1 ± 1.4 t-C ha −1 year −1 and 14.4 ± 1.4 t-C ha −1 year −1 in 2021. The reduction in the net CO 2 emission occurred because GPP (from zero to −12.1 ± 0.7 t-C ha −1 year −1 ) increased more than R eco (from 6.6 ± 1.4 t-C ha −1 year −1 to 13.8 ± 1.1 t C ha −1 year −1 ) between 2018 to 2021 in both areas. The interannual variations in the modelled annual CO 2 flux components were similar between the control and DC treatment areas over the four years. Averaged over the four study years, there was no significant difference between the model estimates of annual NEE in the control and DC areas.  The interannual variations in the modelled annual CH4 balance over the four years differed between the control and DC treatment areas. Specifically, in the control area, the annual CH4 flux increased from −2.2 to 27 kg C ha −1 year −1 between 2018 and 2021. Within the control area, the observed CH4 emissions at 4 m from the uncleaned ditches were by 1 and 2 magnitudes higher relative to the 20 and 40 m locations, respectively, leading to a contribution of 6% and 25% to the carbon balance at this location in 2020 and 2021, respectively. In the DC area, however, the inter-annual variation in the CH4 balance was relatively small and stable within a range of −3.0 to −2.3 kg-C ha −1 year −1 , which represented < The interannual variations in the modelled annual CH 4 balance over the four years differed between the control and DC treatment areas. Specifically, in the control area, the annual CH 4 flux increased from −2.2 to 27 kg C ha −1 year −1 between 2018 and 2021. Within the control area, the observed CH 4 emissions at 4 m from the uncleaned ditches were by 1 and 2 magnitudes higher relative to the 20 and 40 m locations, respectively, leading to a contribution of 6% and 25% to the carbon balance at this location in 2020 and 2021, respectively. In the DC area, however, the inter-annual variation in the CH 4 balance was relatively small and stable within a range of −3.0 to −2.3 kg-C ha −1 year −1 , which represented < 0.15% of the total carbon balances in all four study years. However, when considering the global warming potential of CH 4 as 86 times higher relative to CO 2 over a 20-year time frame (IPCC, 2014), the median contribution of CH 4 to the annual GHG budget (in t CO 2 -eq ha −1 year −1 ) increased from 0.1% in 2018 to 39% in 2021 in the control area, but remained within <5% in the DC area, respectively. The enhanced CH 4 emissions resulted in a 23% increase in the total annual GHG emissions in the control area in 2021; however, given the large uncertainty (from the spatial variation) associated with the annual CH 4 flux estimate, the total annual GHG balances were not significantly different between the control and DC areas.

Ditch Cleaning Effects on Hydrology, Vegetation and GHG Fluxes
This study investigated DC's effects on CO 2 and CH 4 fluxes in a two-way experimental set up, i.e., through (i) a direct comparison of fluxes in the control and DC area, respectively, and (ii) a pre-versus post-treatment assessment in the DC area. Both evaluation alternatives indicated only a minor effect of DC on the net CO 2 exchange and its component fluxes GPP and R eco. Instead, both spatial and temporal variations were mainly dependent on the within-site variations in herbaceous ground vegetation during the post-harvest years. The similar spatial pattern of herbaceous ground vegetation before and after DC indicated the minor impact of DC on the ground vegetation development. In addition, GPP and R eco were strongly correlated with variations in soil carbon and nitrogen, which is in line with previous studies suggesting that soil fertility and organic content are the primary controls over spatial variations of CO 2 fluxes [49,50]. This is because the decomposition of organic matter and associated nutrient release are tightly coupled to microbial respiration and plant carbon uptake [51]. Essentially, positive feedback may develop as the increased plant growth provides additional organic matter input, thereby speeding up the carbon cycle [50]. The spatial variations in both ground vegetation together with soil carbon and nitrogen at our study site itself did not correspond to DC effects and instead might be a relict of the historical drainage effects and/or disturbance during harvest and site preparation.
The lack of WTL changes in the DC area after DC appeared surprising at first. However, this could be explained by the increase in precipitation during the post-DC years, which may have counterbalanced the enhanced drainage function following DC. In fact, the concurrent rise in WTL (by 14 cm) in the control area reveals an increased drainage function in the DC area during the two post-DC years. Given the wetter conditions, the control area was characterised by a relatively shallow WTL (−17 cm on average) during these two post-harvest years, resulting in a narrow upper oxic zone which likely supressed decomposition and the growth of vascular plants, thus explaining the limited R eco and GPP noted in the control area during these wet years (Table 2). However, we observed no significant correlation between WTL and CO 2 flux components (i.e., R eco and GPP) during periods when WTL was <−20 cm, which agrees with similar findings from a peatland forest clear-cut in boreal Finland [26]. This indicates that DC's effects on the Reco, GPP and NEE largely depend on both the initial WTL and the effectiveness of the DC, which jointly regulate the change in soil surface oxic conditions and subsequently the production and decomposition of organic content in the initial years following harvest.
Contrary to our expectations, we did not observe a clear pattern in the effect of ditch distance on WTL, soil moisture and CO 2 fluxes. Given the constantly shallow WTL condi-tions (−17 cm on average during the two post-harvest years) at all distances from uncleaned ditches, the R eco and GPP were higher at locations closer to these uncleaned ditches, which could be due to the contribution from the more extensive growth of Sphagnum moss under the wet conditions. In the DC area, however, the highest R eco and GPP values observed at 20 m from the cleaned ditches may have resulted from optimum WTLs of −35 to −40 cm that created a favourable balance between sufficient water and oxygen availability for enhancing vascular plant growth and decomposition processes [52,53]. Such an impact could have dominated over the drainage legacy effects, as observed in the control area with limited WTL variations. Since the year-to-year increase in GPP was greater than that of R eco at 20 m from the cleaned ditches, the most negative measured daytime NEE (−156 ± 39 mg C m −2 h −1 ) indicates that such a WTL condition likely favours net CO 2 uptake. Nevertheless, we cannot rule out the possibility that the high rainfall during both post-DC years may have influenced the moisture condition and eventually the locations where the optimal WTL with highest fluxes were observed. Thus, additional observations during dry years will be required to further elucidate the interactions among ditch distance, soil water levels and CO 2 fluxes. CH 4 fluxes are known to be strongly related to soil oxic conditions [54,55] and thus are likely modified by a reduction in WTL following DC. Specifically, the deeper mean WTL of −35 ± 2 cm in the DC area retained a substantial surface oxic layer and thereby provided favourable conditions for CH 4 oxidation while limiting CH 4 production to the deeper soil layer [26]. In comparison, the majority (75%) of plots in the control area, including all plots at 4 m distance from the ditches, were characterised by relatively high WTL (between 0 and −20 cm) and close to saturated soil conditions (> 0.45 m 3 m −3 ), which may have facilitated large CH 4 production and emissions from the enhanced anaerobic soil layer [26,56,57]. However, it is noteworthy that CH 4 emission spikes occurred more frequently at several DC plots given the same or even deeper WTL during the two post-DC years. This is likely due to the enhanced amount of vascular ground vegetation (e.g., Deschampsia flex spp.) noted in these DC plots, with roots reaching deep into the CH 4 production zone, providing substrate from root exudates to methanogens and supporting the plant-mediated transport of CH 4 into the atmosphere [58][59][60]. Alternatively, more frequent rain events might have resulted in additional CH 4 production from waterlogged microsites within the oxic soil layer [61]. Altogether, these findings highlight that, relative to the CO 2 exchange, CH 4 fluxes were more sensitive to DC via its effects on soil hydrological properties and plant growth.

Ditch Cleaning Effects on the Annual C and GHG Balances of a Boreal Forest Clear-Cut
We observed a rapid decrease in annual carbon emissions during the four measurement years, which indicates that the high emissions commonly expected to occur in boreal clear-cuts [26,62] might be limited to only the initial few years. One reason for this decrease was the rapid development of the herbaceous ground vegetation, which resulted in increasing GPP. The importance of the fast Recovery of ground vegetation was previously highlighted in nutrient-poor boreal forest clear-cuts on mineral soils [63][64][65]. In comparison, the increase in R eco over the four years was relatively smaller, possibly because the wet conditions in 2020/21 might have constrained microbial aerobic decomposition. Previous studies have shown that the Recovery rate of the net carbon balance after clear-cutting depends on regeneration management, biomass production and decomposition rate, suggesting that it could take anywhere between 1 and 20 years for GPP to counterbalance R eco [62,[66][67][68][69].
Despite the higher daytime R eco and GPP measured in the DC area, our model estimates of their annual sums were not different between the two experimental areas during both post-DC years. The apparent discrepancy in the treatment effect between the modelled annual sums relative to the instantaneous daytime measurements is explained by the greater vegetation coverage in the surrounding area compared to the measurement frames in the control treatment (Figure 4). While the measured fluxes corresponded to the vegetation coverage within the frames, the vegetation coverage in the surrounding area was used as input in the model extrapolation with the goal to obtain estimates that were representative of the entire area. The higher vegetation coverage in the surrounding control area relative to the DC area was primarily a result of the higher prevalence of Sphagnum moss species under favourable shallow WTL conditions. We note that our results are confined to the initial vegetation responses after DC, whereas continued vascular species and seedling development will likely result in a greater contribution from plant carbon uptake to the total respiration at the DC area in the future.
Our result that CH 4 had a minor (<2.5%) contribution to the total carbon balance is in agreement with several previous studies conducted in other boreal forest clear-cuts on drained mineral soils [62,63,70] as well as with findings from peatland forest clear-cuts [26]. However, in terms of the total GHG balance where, the GWP of CH 4 is 86 times over CO 2 , CH 4 emissions become increasingly important when the CO 2 balance is closer towards carbon neutral, i.e., during the later phase when the site shifts from an annual source to sink (Figure 7). Leaving the ditches uncleaned might further create local hotspots for CH 4 emission within the waterlogged areas near the ditch. In comparison, DC effectively mitigated CH 4 emissions by maintaining low WTLs deeper than −30 cm during the wet post-harvest years, resulting in even slightly negative (i.e., uptake) annual GHG balance. This is in agreement with a previous study suggesting that post-harvest CH 4 emissions increase only if the WTL rises above approximately −20 cm [65]. It is noteworthy that the modelled annual CH 4 flux sums were characterised by considerable uncertainty due to the high spatial variability inherent to our measurement plots. The use of the eddy covariance technique is required to overcome this challenge and to provide more precise estimates for the CH 4 balance of heterogeneous clear-cut areas [71].
We caution that our study lacked flux data from the snow-covered winter periods, and our annual flux budget therefore relied on simple assumptions for estimating the contribution from winter fluxes. However, previous studies based on continuous year-round data collected with the eddy covariance and automated chamber techniques suggest low fluxes during snow-covered winter periods commonly contributing to less than 10% of the annual CO 2 and CH4 budgets in boreal ecosystems [62,72,73]. This indicates a limited potential bias introduced by our assumptions in estimating the winter flux. It is further noteworthy that this study only addressed the initial responses of the environment and fluxes to DC. However, long-term effects on the growth of the tree layer might further enhance the DC impact on the forest carbon balance over an entire rotation period. Moreover, the ditch functions will again deteriorate over time through erosion and sedimentation, which may result in a nonlinear trajectory of DC effects on the soil water dynamics and thus carbon balance. Thus, more empirical data from winter periods and extending over the decadal timescale of post-DC years are required to better understand the DC effects on the forest ecosystem carbon balance across contrasting sites. Given the steady increase in DC activities, particularly within Fennoscandia, it is essential to improve our knowledge of DC effects on the forest ecosystem-atmosphere exchange of carbon and GHGs with the goal to provide an empirical knowledge base that can support the development of sustainable and climate-responsible forest management strategies.

Conclusions
This study investigated how the cleaning of degraded drainage ditches affects CO 2 and CH 4 fluxes from a drained forest clear-cut on mineral soil in boreal Sweden over four years (i.e., 2-5 years after harvest). Based on our findings, we conclude that: (1) The clear-cut area with old and degraded ditches acted as a net carbon source in all four post-harvest years. However, during the study period, the annual total carbon emissions decreased by 76% (from 6.7 ± 1.4 t-C ha −1 year −1 to 1.6 ± 1.6 t-C ha −1 year −1 ).
(2) Ditch cleaning had a limited initial effect on the spatio-temporal variations in the net CO 2 exchange and its component fluxes, GPP and R eco. The variation in the component fluxes was instead primarily controlled by the within-site variations in ground vegetation development, likely in response to drainage legacy effects on soil carbon and nitrogen contents. (3) In comparison, ditch cleaning reduced the soil water content and thereby mitigated CH 4 emissions during wet post-harvest years. (4) Ditch distance had no consistent effect on CO 2 and CH 4 fluxes. While R eco and GPP tended to increase towards uncleaned ditches coinciding with legacy trends in soil carbon and nitrogen content, maximum R eco and GPP occurred at 20 m from cleaned ditches, likely in response to an optimal WTL for vascular plants commonly observed in the DC area. In comparison, high emissions of CH 4 mainly occurred on nearly saturated soil locations near to uncleaned ditches. (5) Overall, ditch cleaning had no significant impact on the annual carbon and GHG balance in the initial post-harvest years. (6) There is a critical need for long-term observations to evaluate DC effects on the forest carbon and GHG balances during the entire subsequent rotation period.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/f13060842/s1, Table S1: Annual mean ± standard error (SE) and mixed effect model results for environment variables during manual flux measurements during pre-DC (2018/19) and post-DC periods (2020/21), classified by ditch treatment effects (Control versus DC) and distance to ditches combinations (4m, 20m and 40m). Fixed factors of mixed effect models include clean treatment (T), distance to ditches (D) and their interaction (TD). Column N refers to the sample size of each model. Significant p; Table S2: Model estimates for the total annual carbon dioxide (CO 2 ) and methane (CH 4 ) balances in control and ditch cleaning (DC) treatment areas. CO 2 fluxes includes its net exchange (NEE) and component fluxes of gross primary productivity (GPP) and ecosystem respiration (R eco ). Columns denotes the three distances (4m, 20m and 40m) from ditches estimated for the four study years (2018 to 2021). Values are in the unit of t-C ha −1 year −1 for CO 2 and of kg-C ha −1 year −1 for CH 4. Numbers are represented with ± standard error (SE); Table S3: The correlation coefficients (or loadings) of the each significant Principal components (PCs) with the input variables. Significant PCs are defined by the Kaiser criterion where the eigenvalues are greater than one. Five PCs were defined significant in the PCA with sample data for Pre-DC Period shown at the left columns, whereas four PCs were significant for the Post-DC period at the right columns. The total variance explained for each PC (in %) is presented at the first row. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Measured data used in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Additional Information on the Measurement of N 2 O Fluxes
N 2 O fluxes were measured biweekly with a separate set of opaque chambers (48.5 × 48.5 × 50 cm) in 2020. The chambers were placed on the frames for 75 min, during which, four 60 mL gas samples were taken from the chamber with plastic syringes at 0, 25, 50 and 75 min after closure. During chamber closure, a small fan was operated inside the chamber to maintain air circulation, and a Hobo ® pendant temperature logger (Onset Computers, Bourne, MA, USA) was provided to continuously Record air temperature at 5 s intervals in the chamber headspace. The headspace gas in the 60 mL syringe was then injected into 20 mL evacuated glass vials and determined for their N 2 O concentration within seven days using a headspace sampler (TurboMatrix 110; Perkin-Elmer, Waltham, MA, USA) and gas chromatograph (Clarus 580, PerkinElmer Inc., Waltham, MA, USA). N 2 O was separated by two identical 30 m × 0.53 mm internal diameter megabore capillary porous-layer open tubular columns (Elite PLOT Q) maintained at 30 • C (detection limit: N 2 O < 1 ppb). The GC system was equipped with an electron capture detector (ECD) operated at 375 • C for N 2 O analysis. The linear increase in N 2 O concentrations inside the chamber over time was then converted into a flux estimate using the ideal gas law (Equation (A1)): where F is the measured flux (µmol m −2 s −1 ), dC/dt is the linear slope with the highest r 2 of concentration change over time (ppm s −1 ), V is chamber headspace volume (m 3 ), p is the atmospheric pressure (approximated by a constant value of 101,325 Pa), R is the universal gas constant of 8.3143 (m 3 Pa K −1 mol −1 ), T a is the mean air temperature (K) during the measurement, and A is the frame area (m 2 ).