Initial effects of post-harvest ditch cleaning on greenhouse gas fluxes in a hemiboreal peatland forest

Ditch cleaning (DC) is a well-established forestry practice across Fennoscandia to lower water table levels (WTL) and thereby facilitate the establishment of tree seedlings following clear-cutting. However, the implications from these activities for ecosystem-atmosphere greenhouse gas (GHG) exchanges are poorly understood at present. We assessed the initial DC effects on the GHG fluxes in a forest clear-cut on a drained fertile peatland in hemiboreal Sweden, by comparing chamber measurements of carbon dioxide (CO 2 ), methane (CH 4 ) and nitrous oxide (N 2 O) fluxes from soil and ditches in DC and uncleaned (UC) areas over the first two post-harvest years. We also evaluated spatial effects by comparing fluxes at 4 m and 40 m from ditches. We found that 2 years after DC, mean ( ± standard error) WTL of (cid:0) 65 ± 2 cm was significantly lower in the DC area compared to (cid:0) 56 ± 2 cm in the UC area. We further observed lower gross primary production and ecosystem respiration in the first year after DC which coincided with delayed development of herbaceous ground vegetation. We also found higher CH 4 uptake but no difference in N 2 O fluxes after DC. Greater CH 4 uptake occurred at 4 m compared to 40 m away from both cleaned and uncleaned ditches. Model extrapolation suggests that total annual GHG emissions in the second year were reduced from 49.4 ± 17.0 t-CO 2 -eq-ha (cid:0) 1 -year (cid:0) 1 in the UC area to 27.8 ± 10.3 t-CO 2 -eq-ha (cid:0) 1 -year (cid:0) 1 in the DC area. A flux partitioning approach suggested that this was likely caused by decreased heterotrophic respiration, possibly because of enhanced soil dryness following DC during the dry meteorological conditions. CH 4 and N 2 O fluxes from clear-cut areas contributed < 2 % to the total (soil, ditches) GHG budget. Similarly the area-weighted contributions by CO 2 and CH 4 emissions from both cleaned and uncleaned ditches were < 2 %. Thus, our study highlights that DC may considerably alter the post-harvest GHG fluxes of drained peatland forests. However, long-term observations under various site conditions and forest rotation stages are warranted to better understand DC effects on the forest GHG balance.


Introduction
Northern peatlands represent an important sink for atmospheric carbon dioxide (CO 2 ) and constitute a major global store of soil carbon (270-621 Gt C) (Adams andFaure, 1998, Loisel et al., 2014;Turunen et al., 2002;Yu, 2012). However, during the past century large areas of pristine mires were drained for forestry activities in Fennoscandia, the Baltics, Canada, USA and Russia (Minkkinen et al., 2008;Paavilainen and Päivänen, 1995;Päivänen and Hånell, 2012). Today, many forests on drained peatlands have reached the end of their rotation age and are subject to harvest. Following final harvest (i.e. clear-cutting), ditch cleaning (DC) is a common practise in northern Europe (Finland, Sweden, Estonia and Latvia) to improve drainage and thereby facilitate tree seedling establishment (Bergquist et al., 2016;Finér et al., 2018;Päivänen and Hånell, 2012;Tomppo, 2005). According to The Swedish National Forest Inventory (NFI), the 5-year average rate of forest area subject to DC activities in Sweden has increased by a factor of four, from 2600 ha year − 1 in [2005][2006][2007][2008][2009] to 10400 ha year − 1 in 2015-2019 ( Fig. A.1). From 2003 to 2019, ~102 000 ha have been ditch cleaned in Sweden, representing about 10 % of the drained productive forest area (Hånell and Magnusson, 2005). The increasing trend of DC activities in recent years is particularly pronounced for forests on formerly drained peatlands.
While DC may be a beneficial measure for supporting seedling establishment and tree growth, a deeper water table level (WTL) after DC may strongly affect soil biogeochemistry and vegetation growth, resulting in multiple and complex consequences for greenhouse gas (GHG) production and consumption processes. Specifically, drainage on wet soils may accelerate peat decomposition rates due to increased oxygen availability (Drzymulska, 2016) leading to higher soil carbon dioxide (CO 2 ) emissions (Maljanen et al., 2010;Ojanen et al., 2013;van Huissteden et al., 2006). At the same time, drainage of wet soils increases root aeration and nutrient availability, which facilitates establishment of initial ground vegetation as well as tree seedlings and subsequent tree growth, enhancing rates of gross primary productivity (GPP) and carbon sequestration (Hökkä and Kojola, 2001;Hökkä and Kojola, 2003;Lauhanen and Ahti, 2001;Sikström et al., 2020;Sikström and Hökkä, 2016).
Following drainage, methane (CH 4 ) emissions can be expected to decrease due to improved soil aeration which may potentially turn the soil into a net CH 4 sink (Kasimir et al., 2018;Maljanen et al., 2001;Martikainen et al., 1995;Nykänen et al., 1998;von Arnold et al., 2005). Yet the extent of such decrease and the thresholds regulating the switch in the source-sink function remain unclear, due to partly compensating effects. Specifically, on the one hand, drainage extends the oxic zone and increases the potential for aerobic CH 4 oxidation while forcing the zone for CH 4 production deeper into the soil (Borken et al., 2006;Feng et al., 2020;Fest et al., 2017). At deeper soil profiles, methanogenesis could be further limited by substrate supply as decomposition rates usually decline with increasing peat age (Frolking et al., 2001). On the other hand, enhanced establishment of vascular plant field layer after DC could increase CH 4 production and emission through enhancing rhizosphere exudation of substrates and CH 4 transport via aerenchyma tissue and stomatal conductance (Chu et al., 2014;Garnet et al., 2005;Granberg et al., 1997;Long et al., 2010). Interaction effects among these separate controls, e.g., whether rooting depth extends below the mean WTL might further complicate the net effect of DC effects on the net CH4 exchange.
Drained organic soils are often characterized by low C:N ratios (Ernfors et al., 2007). Under high soil nitrogen availability, i.e., C:N ratio < 20, considerable production and emission of nitrous oxide (N 2 O) may occur through interacting biological pathways for reducing NH 4 + and NO 3 − (Firestone and Davidson, 1989;Klemedtsson et al., 2005). WTL drawdown is expected to interfere such pathways through, for instance, triggering N 2 O emissions through incomplete denitrification under partially-oxidised conditions (Rubol et al., 2012) and increasing nitrifier activity when the soil is further aerated (Santin et al., 2017). However, further drainage under dry conditions could possibly reduce N 2 O emission by limiting denitrifying bacterial activities (Christiansen et al., 2012;Rassamee et al., 2011). Thus, N 2 O production responds to soil moisture along an optimum-curve which explains both positive (e.g. Pihlatie et al., 2004;Rochette et al., 2010) and negative (e.g. Christiansen et al., 2012;Pärn et al., 2018) correlations between N 2 O production and WTL in previous studies. Thus, the net effect of DC on N 2 O emissions likely depends on the combination of the initial WTL, the effectiveness of the DC measure in lowering WTL, and soil nutrient status. At present, empirical data exploring these complex relationships are lacking. Given the potential increase in N 2 O emissions following clear-cutting of drained peatland forests (Huttunen et al., 2003), it is critical to improve our understanding on the response of N 2 O emission to management effects such as post-harvest DC. Drainage ditches are potentially important GHG emission hotpots and contributors to ecosystem-scale GHG budgets. Ditches can be significant CH 4 sources due to their commonly anoxic soil conditions which stimulate methanogenesis (Hyvönen et al., 2013;Minkkinen and Laine, 2006;Peacock et al., 2017;Sundh et al., 2000). Ditches also exchange CO 2 (e.g. Sundh et al., 2000;Teh et al., 2011;Hyvönen et al., 2013;Vermaat et al., 2011) and may emit N 2 O (e.g. Reay et al., 2003;Teh et al., 2011;Hyvönen et al., 2013). DC could reduce the retention of suspended solids in ditches , which consequently may decrease the availability and quality of substrates for facilitating GHG emissions (Hyvönen et al., 2013). To date, however, studies of CO 2 and N 2 O fluxes from drainage ditches are too few to draw reliable conclusions (Evans et al., 2016). It is thus important to improve our understanding of DC effects on ditch GHG fluxes and to evaluate their contribution to ecosystem-scale GHG budgets.
Earlier studies of effects of natural WTL variations on forest GHG fluxes (e.g. Korkiakoski et al., 2017;Pearson et al., 2012) indicate that DC might alter the GHG balance significantly through its impact on ecosystem hydrology. While DC effects on the GHG balance have been recently explored for a forest clearcut on wet mineral soil in boreal Sweden (Tong et al., 2022), such assessment is currently lacking for the nutrient-rich peat soils in the hemiboreal region. Based on a paired experimental setup including adjacent DC and UC areas, the aim of this study was therefore to investigate the impact of post-harvest DC on GHG fluxes from clear-cut areas and ditches in a drained fertile peatland in hemiboreal Sweden. The specific objectives were to: (1) investigate initial effects of DC following clear-cutting on spatio-temporal variations in GHG fluxes; (2) identify environmental factors that drive changes in GHG fluxes in response to DC; and (3) assess DC effects on the annual GHG balances including fluxes from clear-cut areas and ditches.

Site description and experimental design
This study was conducted at the Skogforsk experimental field site 303 Tobo which is located approximately 40 km north of Uppsala city, Sweden (60 • 15 ′ 54 ′′ N, 17 • 37 ′ 00 ′′ E, 40 m.a.s.l.). Climate normals for 1991-2020 recorded at a nearby (~15 km) weather station (Film) by the Swedish Meteorological and Hydrological Institute (SMHI; https://www .smhi.se), had annual mean air temperature and precipitation of 6.2 • C and 603 mm respectively, and a snow cover period (mean snow depth at least 5 cm) lasting from November to April. The study area was originally a minerotrophic mire which was drained with a ditch network for agricultural use (Fig. 1). In the beginning of the 1970 s Norway spruce (Picea abies) seedlings were planted at the site. According to the national SLU Soil Moisture Map (Ågren et al., 2021), the hydrological site conditions at Tobo in 2011 (i.e., prior to harvest) were similar to those of other nearby peatland forest areas (Fig. A.2).
In September 2017, the ca. 45-year-old stand consisting of about 90 % Norway spruce and 10 % birch (Betula pendula) was harvested. In July 2016, a total of 11 main experimental plots were installed as a splitplot experiment with ditch cleaning as the main treatment (UC = uncleaned; DC = ditch cleaned) and site preparation (mounding) as sub-treatment (No mounding or mounding). Each main plot (60 m × 80 m) included two different sub-treatment plots (30 m × 80 m; hereafter subplots) located perpendicular to the ditches (short side parallel to the ditch). In this study, 10 of these 11 subplots with no mounding as sub-treatment were used, with five of them located next to uncleaned ditches (UC ditches) and cleaned ditches (DC ditches), respectively, which were evenly distributed over the area (Fig. 1). Ditch cleaning was undertaken using an excavator during the first week of November 2017. Norway spruce seedlings (15-25 cm in height) were subsequently planted (with a spacing of 2.15 m × 2.15 m) in the beginning of May 2018. Due to high seedling mortality and browsing damage during 2018, the plots were fully replanted in the beginning of May 2019 with similar Norway spruce seedlings and at the same spacing as in 2018. However, replanting of the subplots was carried out only within the first 45 m from the ditch. The ditch network system at the site consists of a central main ditch running north-south across the entire site, in connection with perpendicular tributary ditches, four of which are situated east of the central main ditch at about 185 m intervals, and two situated west of the central ditch. There is also a drainage ditch at the southern edge of the site (Fig. 1). DC was performed in the central ditch, two of the eastern perpendicular tributary ditches and both western tributary ditches (Fig. 1). The site lies on a relatively flat area, but has a gentle slope (~0.3 %) due to which drainage water flows southward along the central ditch.
The field-vegetation type according to Hånell (1986), assessed in August 2019, was predominantly classified as a tall-herb type indicating a site quality that corresponds to a predicted stem wood production of 10.2 m 3 ha − 1 year − 1 . Herbaceous species including Chelidonium majus, Cirsium spp., Deschampsia cespitosa, Lactuca muralis, Ranunculus repens, Senecio sylvaticus, Taraxacum sp., Urtica dioica and Vaccinium spp. had established across the site during the first two years following harvest, contributing to a total dry weight of 30 ± 2 ton ha − 1 collected at all 10 subplots in September 2020. Eriophorum vaginatum and Juncus spp. were occasionally found in the uncleaned ditches and likely existed also in the cleaned ditches before DC. The degree of humification of the peat was H8 -H9 (von Post scale;von Post, 1922) at the peat depth of 0 -30 cm. Peat depth was > 100 cm within all subplots used in this study. The mean (±standard error) soil carbon and nitrogen contents determined for all 10 subplots in September 2020 were 50.4 ± 0.6 and 3.06 ± 0.04 % (i.e. C:N ratio: 16.5 ± 0.3), respectively, without significant difference between UC and DC areas (Table A.1 and Fig. A.3). Soil (0-30 cm depth) concentrations of phosphorus, calcium, potassium and magnesium were 0.87 ± 0.03, 53 ± 0.97, 0.17 ± 0.04 and 0.86 ± 0.04 g kg − 1 total solids, respectively, without significant difference between UC and DC areas (Table A.1 and Fig. A.3). While site-specific information is lacking, data from the Swedish Forest Soil Inventory (SFSI) suggest that soil pH ranges between 5.8 and 6.1 in this region. The pH of ditch water samples was measured to be 7.0 ± 0.2. Prior to ditch cleaning (i.e. between the harvest in September to early November 2017), the mean WTL was − 35 ± 3 cm in the UC area and − 33 ± 4 cm in the DC area, revealing statistically insignificant difference (p = 0.74) with each other (Fig. A.3). Thus, although the relatively short period between harvest and ditch cleaning did not allow for extensive pre-treatment measurements, the similar soil chemistry and pre-treatment WTL data across the study site suggest that no bias in local environmental conditions was present that could confound the assessment of ditch cleaning effects.

Greenhouse gas flux measurements
Within each subplot, we selected GHG flux sampling locations at 4 m and 40 m on both sides of the ditches (hereafter fluxes in the clear-cut area) (Fig. 1). In addition, CO 2 and CH 4 were measured in the UC and DC ditches (hereafter fluxes in the ditches). In the clear-cut area, we carried out CO 2 and CH 4 flux measurements every second week using the closed dynamic chamber method (Livingston and Hutchinson, 1995) during daytime from May to November in 2018 and 2019, i.e., during the first two years after clear-cutting and DC. In late April 2018, i.e. one month before the first flux measurement, square aluminium frames (48.5 × 48.5 cm) were permanently installed at each sampling location. The frame base was inserted down to 5 cm below the soil surface. At each frame, measurements were conducted first with a transparent chamber (light transmission rate 81 %) to estimate the net ecosystem CO 2 exchange (NEE). The chamber was then covered with an opaque and light-reflective shroud to estimate the ecosystem respiration (R eco ) during dark conditions. The gross primary productivity (GPP) was then calculated by subtracting R eco from NEE. In the first year (2018), we used a chamber with the dimensions of 48.5 × 48.5 × 30 cm whereas in the second year (2019), a taller chamber (48.5 × 48.5 × 50 cm) of the same material was required to cover the herbaceous vegetation growing within the frames. While the headspace air in the smaller chamber was mixed by the continuous sample air return flow from the analyser (see below), this larger chamber was equipped with a small fan to further support air mixing in the larger and densely vegetated headspace. A comparison of repeated (i.e., within < 5 min) fluxes measured with the two different chamber sizes indicated good agreement for CO 2 (r 2 = 0.95) and CH 4 (r 2 = 0.58) fluxes respectively (Fig. A.4). The lower r 2 for the CH 4 flux was likely caused by the disturbance from the first measurement which may have altered the relatively small CH 4 concentration gradient prior to the subsequent measurements. In any case, while the different chamber type might have potentially affected the between-year comparison (i.e., 2018 vs 2019), it had no effect on the assessment of DC effects since the same chamber was used in both treatments within each of the two years. During each measurement, the chamber was placed on the frame for 180-240 s. The chamber was connected to an Ultraportable Los Gatos Research (LGR) greenhouse gas analyser (Model 915-0011; San Jose, CA, USA) with a built-in sampling pump drawing and returning air between chamber and analyser at a flow rate of 0.8 standard L min − 1 . The analyser determined chamber CO 2 and CH 4 concentrations at 2 s intervals with accuracy levels of ± 300 ppb and ± 2 ppb, respectively. To eliminate the effect of water vapour on volumetric dilution and spectroscopic line broadening, the dry air mole fraction was used.
The slope of the change in CO 2 and CH 4 concentration over time (dC/dt; ppm s − 1 ) was estimated by a simple linear regression over a chosen data range. Specifically, after discarding the first 20 s concentration data (dead bands), the slopes for all possible 100 s windows (or 60 s if H 2 O has reached saturation during measurement) over the measurement period were calculated and the slope with the highest coefficient of determination (r 2 ) was chosen as dC/dt. The flux rate was then calculated based on dC/dt and the ideal gas law (Eq. (1)): 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 ). N 2 O fluxes were measured every second week with a separate set of opaque chambers (48.5 × 48.5 × 50 cm) from May to November 2019. Each of these chambers was equipped with a small fan for maintaining air circulation and a Hobo® pendant temperature logger (Onset Computers, Bourne, MA, USA) for monitoring air temperature in the chamber headspace. 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. The gas samples were injected into 20 ml evacuated glass vials and analysed for their N 2 O concentration within seven days using a headspace sampler (Turbo-Matrix 110; Perkin-Elmer, MA, USA) and a gas chromatograph (GC) system (Clarus 580, PerkinElmer Inc, USA) fitted with two identical 30 m × 0.53 mm inners 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 of N 2 O concentrations inside the chamber over time was then converted into a flux estimate using Eq. (1).
Poor quality flux data, defined by the root-mean-square error (RMSE) and r 2 of the chosen slope of dC/dt, were filtered out before further analysis. Specifically, in the clear-cut area, CO 2 fluxes with RMSE > 2.5 ppm and r 2 < 0.90, CH 4 fluxes with RMSE > 2.5 ppb and r 2 < 0.90, and N 2 O fluxes with RMSE > 10 ppb and r 2 < 0.90 were removed. These thresholds were identified based on visual examination of the data (Fig. A.5). These quality control procedures led to the removal of about 2 %, 5 %, and 8 % of all CO 2 , CH 4 and N 2 O fluxes measured, respectively. Note that the sign convention in this study is such that positive and negative values indicate that the ecosystem is a source and sink, respectively.
Ditch CO 2 and CH 4 fluxes were measured monthly with opaque floating chambers (diameter 31.5 cm, volume 9.56 L) in 2018 and 2019 from seven ditches, using a Picarro G4301 GasScouter (Picarro, Santa Clara, CA, USA) at a sampling interval of 1 s and accuracies of ± 400 ppb and ± 3 ppb for CO 2 and CH 4 concentrations, respectively. Ditch water levels were recorded simultaneously during each measurement using a ruler. On some occasions when the ditches were dry, the chamber was gently pushed into the sediment to create a seal as commonly applied in recent studies (Peacock et al., 2021a;Peacock et al., 2021b). In addition, since it was not possible to fit the floating chamber over some of the tall vegetation species (e.g., Juncus spp.) growing sporadically in some ditches, we were not able to account for their GPP and autotrophic respiration (R a ). Thus, our ditch measurements represent the ditch water (or sediment) surface-atmosphere flux. Low quality ditch fluxes were removed with the threshold RMSE > 6 ppm and r 2 < 0.90 for CO 2 fluxes, and RMSE > 5 ppb and r 2 < 0.90 for CH 4 fluxes ( Fig. A.6). Ditch N 2 O fluxes were not determined in this study, however, their magnitudes are commonly very low (Peacock et al., 2017).

Environmental variables
Environmental conditions were recorded both manually during each flux sampling campaign and continuously in hourly intervals through the entire year. The manually recorded data was used to i) investigate the environmental drivers of the measured fluxes through statistical analyses and ii) to calibrate the continuous data which served as input for the models to estimate annual fluxes (See section 3.5). Manual WTL measurements were taken inside PVC groundwater tubes (Ø = 32 mm external and 26 mm inside, 125 cm long with 3 mm holes every 2.5 cm) adjacent to each measurement frame and inserted to about 1 m depth into the peat. Air temperature (T a ) along with soil temperature at 5 and 10 cm depth (T s5 , T s10 ) outside the frame were recorded manually during each flux measurement using a handheld temperature meter (shaded from direct sunlight during T a measurement). In 2019, photosynthetically active radiation (PAR) was measured with a handheld radiometer (QSO-S PAR Photon Flux Sensor connected with ProCheck data logger, both by Decagon Devices, Pullman, WA, USA) and soil moisture within the upper 5 cm (SM) was measured at three sides around the frame during each flux measurement using a GS3 combined moisture-temperature sensor (Decagon Devices, Pullman, WA, USA) connected to the ProCheck data logger.
Continuous hourly data as input for the models to estimate annual GHG fluxes were obtained at the site and from nearby weather stations. Specifically, the WTL was monitored at an hourly interval using WT-HR 1000 probes (TruTrack ltd, Christchurch, New Zealand) placed inside the PVC groundwater tubes adjacent to each frame. In addition, automated soil temperature and moisture sensors (CS655; Campbell Scientific, Logan, USA) connected to a data logger (CR1000X; Campbell Scientific, Logan, USA), were installed at one location near the centre of the study area to monitor the temporal variations at hourly intervals ( Fig. 1). Continuous hourly PAR and T a data were obtained from the nearest available weather stations ICOS-Norunda (~21 km away) and the SMHI meteorological station in Film (~15 km away), respectively. For all environmental variables (i.e., WTL, PAR and T a ), the linear correlations between manual and automated sensor data were strong (r 2 > 0.90, Fig. A.7). The continuous environmental data were therefore calibrated with the manual measurements to obtain hourly model input data adjusted to the study site conditions.

Vegetation data
To assess the effects of vegetation development on GHG fluxes, we determined the ground vegetation areal coverage defined as the projected area of vegetation over a unit of land (m 2 m − 2 ). In addition, we took overhead images of each frame in October 2018 and in July 2019 to derive a vegetation greenness index defined by the green chromatic coordinate (g cc ) (Järveoja et al., 2016a;Peichl et al., 2015;Sonnentag et al., 2012) (Eq. (2)).
where g cc refers to the greenness index from the image taken on the frame; R, G and B denote intensity (0-255) of the red, green and blue image channels. The RGB values were averaged for each image pixel located within the chamber frame. To describe the growing season phenology as model input for estimating annual GHG balances in 2019, a Gaussian curve was fitted as a function of number of days away from the assumed day of peak vegetation growth Wilson et al., 2007) where g cc (JD) denotes the greenness index on the particular Julian day (JD; day of year numbered from 1 to 365); g ccmax refers to the maximum greenness index (see Eq. (2)) which was estimated to be reached in early July 2019 (JD max = 189 on visual field observations of plant development throughout growing season; Parameter b denotes the width of the curve. In July 2019, vegetation growth outside the frames was also examined by measuring ground vegetation height, areal coverage and g cc at 12 spots evenly distributed within 15 m around each frame, in order to provide information on (1) whether the vegetation inside the chamber frames was representative of the surrounding area, and (2) DC effects on ground vegetation growth.

Modelling of annual GHG budgets
To estimate the annual CO 2 and CH 4 balances, we developed nonlinear regression models to predict hourly fluxes in response to environmental parameters and vegetation development following Järveoja et al. (2016aJärveoja et al. ( , 2016b, Kandel et al. (2013) and Olson et al. (2013). Particularly, GPP from each frame was fitted to hourly mean PAR using a hyperbolic function modified with normalised framespecific g cc which represents seasonal changes in vegetation biomass (Eq. (4)): where GPP denotes the hourly gross primary production (mg m − 2 h − 1 of CO 2 -C); α denotes model fitted value of the initial slope of the light-use efficiency of photosynthesis (mg µmol photons − 1 of CO 2 -C); PAR denotes the hourly mean photosynthetically active radiation (µmol m − 2 s − 1 ); P max denotes the modelled fitted value of maximum photosynthesis under light saturation (mg m − 2 h − 1 of CO 2 -C); and g ccnorm is the daily estimated frame-specific chromatic greenness index (g cc (JD)) normalized to scale between 0 and 1.
In the R eco model, we used an exponential relationship with T a based on Lloyd and Taylor (1994) modified with the normalised frame-specific g cc , as second explanatory variable (Eq. (5)): which denotes sensitivity of respiration to T a , and β which is a scaling parameter for plant development representing the contribution of plant autotrophic respiration (R a ) to R eco . Using continuous hourly mean T a and ambient PAR data from the nearby weather station as well as g cc of the frame surrounding area as input variables to the respective models, diel hourly R eco and GPP were modelled and summed up for the entire year.
To estimate annual autotrophic and heterotrophic respiration (R a and R h ) contributions to R eco , we assumed a carbon use efficiency (CUE) of 0.5 to derive R a = GPP × CUE and R h = R eco − R a (Waring et al., 1998;Gifford, 2003). During winter periods (i.e., November to April) the model estimated low-temperature fluxes mainly based on R 0 . It is further noteworthy that modelling nighttime respiration based on the nighttime temperature assumes a constant diel temperature sensitivity.
Hourly estimates CH 4 fluxes were modelled using an exponential relationship with WTL and soil temperature at 10 cm depth (Olson et al., 2013) (Eq. (6)): where CH 4 denotes hourly CH 4 flux (g m − 2 h − 1 of CH 4 -C); b 1 and b 2 denote the model fitted sensitivity of CH 4 flux to water table level (WTL, cm) and soil temperature at 10 cm depth (T s10 , • C), respectively; b 0 denotes the intercept of the model. Site-level data of continuous WTL and T soil were used as input for modelling hourly CH 4 fluxes. We applied Monte Carlo simulations to evaluate the uncertainty of the annual flux budgets estimated by the model extrapolations (Smith and Heath, 2001). For this purpose, a normal distribution was assigned to each model input parameter based on its mean and standard deviation derived during model development (Table 1). Then, a large number (1000) of possible scenarios were generated using random values from the normal distribution of each model parameter. The standard deviation for the set of 1000 predicted annual sums was then used to describe the uncertainty of the annual flux budget estimates.
Separate models were developed for UC and DC areas with the determination of coefficient (R 2 ) as the criteria for selecting the best final models. The model parameters for the two areas and study years are summarized in Table 1. The comparisons between measured and model fitted values from all nonlinear models are presented in Fig. A.8. Due to the weak response of soil N 2 O as well as ditch CO 2 and CH 4 emissions to environmental variables, annual sums of these fluxes were estimated by scaling the median of the measured fluxes to the entire year. Modelled estimates of CH 4 and N 2 O were transformed into CO 2equivalents (CO 2 eq) by applying their global warming potential (GWP) of 34 and 298 over a 100 year timeframe, respectively (IPCC, 2013). The median rather than the mean of the measured CH 4 and N 2 O fluxes was used to avoid overestimating the annual sum due to episodic high fluxes.

Statistical analysis
We first performed a principal component analysis (PCA) to explore the overall coherence structures among all the variables, including the CO 2 and CH 4 fluxes, environmental factors and treatments (DC and distance). The relationship of N 2 O fluxes with potential factors were analysed in separated PCA analysis. The input variables to the PCA were normalized to zero mean values by subtracting the mean and to unit variance by dividing the values by the standard deviation of the variable (Jolliffe, 1990). Significant principal components were selected and presented using the broken-stick model (Jackson, 1993). The variable loading, defined as the correlation between each variable and PC, was used as a criterion to determine the relationships (Cadima and Jolliffe, 1995).
Next, since PCA did not provide information on the significance level on the treatment effect on the study variables, we applied mixed effect models with repeated measures to quantify the statistical significance level of the treatment effects (i.e., DC treatment and distance to ditch) on the spatio-temporal variation of environmental (i.e., T a, T s10 , T s5 , SM or WTL) and GHG flux (i.e. GPP, R eco , NEE, CH 4 or N 2 O) variables (Eq. (7)).

Table 1
Model parameters for estimating gross primary production (GPP) (Eq. (4)), ecosystem respiration (R eco ) (Eq. (5)) and methane (CH 4 ) fluxes (Eq. (6)) for uncleaned (UC) and ditch cleaned (DC) clear-cut area in the study year 2019; α is the initial slope of the light-use 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 air temperature (T a ); β represents the contribution of autotrophic respiration to R eco ; b 0 denote the intercept of the CH 4 function, b 1 and b 2 are the sensitivity of CH 4 flux to water table level and soil temperature at 10 cm depth, respectively; numbers in parentheses indicate standard error; R 2 denotes the coefficient of determination of the model, RMSE and MAE denotes the rootmean-square error and mean absolute error of the model with the same unit as the predicted fluxes (i.e. mg m − 2 h − 1 for CO 2 and mg m − 2 h − 1 for CH 4 ). These models included a spatial covariance structure where correlations decline over time (Phillips et al., 2001). The statistical models applied were as follows: where y ijk denotes the environmental or GHG flux variable for sampling occasion i with DC treatment j (j = UC or DC) at distance to ditch k (k = 4 m or 40 m); μ denotes the overall mean of the environmental or GHG flux variable; T j denotes the fixed effect of DC treatment j; D k denotes the fixed effect of distance to ditch k; TD jk denotes the two-way interaction between the effects of the treatment j and distance to ditch k; S ijk denotes the random effect of sampling occasion i; ε ijk denotes the random error for sampling occasion i with treatment j at distance to ditch k. Mixed effect models were proven robust to different data distributions (Schielzeth et al., 2020). Statistical results from the mixed effect models were considered significant at p < 0.05. The standard error (±SE) of the sample averages was used as a measure of uncertainty throughout this paper. All statistical analysis was conducted using the Mathworks Matlab software R2019b.

Environmental data
The annual mean air temperatures in 2018 and 2019 at the nearby SMHI weather station were 6.8 and 6.7 • C, which are close to the 30year long-term average (1991-2020) of 6.2 • C. However, May 2018 was an unusually warm month based on both local site data (Fig. 2) and SMHI weather station data (2018: 14.3 • C; 1991-2020: 10.0 • C). The annual precipitation in 2018 (522 mm) and 2019 (703 mm) measured at the SMHI weather station was below and above the long-term average 30-year average (606 mm), respectively. Except during July and August 2018 with 37 % to 41 % higher precipitation amounts than the 30-year normal, the monthly precipitation sums during 2018 were 21 % to 83 % less than the 30-year normal. The temporal dynamics and magnitudes of photosynthetically active radiation (PAR) were similar in both years (Fig. 2).
In 2018, the seasonal variation in the spatial means (i.e., averaging 4 m and 40 m locations) of WTL varied from − 44 to − 62 cm and − 49 to − 65 cm in the UC and DC area, respectively (Fig. 2), with means of − 54 ± 1 and − 57 ± 1 cm (p > 0.05) for the entire monitoring period, May -November (Table A.2). In 2019, the seasonal variation in the spatial mean WTL spanned from − 36 to − 74 cm and from − 42 to − 89 cm in the UC and DC area, respectively, being significantly lower for DC with − 65 ± 2 cm compared to − 56 ± 2 cm for UC (p < 0.01; Table A.2) over the 2019 monitoring period. The lowest WTL occurred from early July to the end of September, during which time the WTL difference between UC and DC sampling locations became largest (varying between 5 and 20 cm). In both years 2018 and 2019, the mean WTL at sampling locations located 4 m from ditches was significantly lower (p < 0.01) than at 40 m distance in both UC and DC areas (Table A.2).
In 2019, similar patterns were observed for soil moisture, with a significantly higher (p < 0.01) growing season mean of 0.34 ± 0.01 m 3 m − 3 in the UC area than 0.30 ± 0.01 m 3 m − 3 in the DC area ( Fig. 2 and Table A.2). The lowest soil moisture level occurred in July and August, when the difference in soil moisture levels between UC and DC sampling locations reached up to 0.11 m 3 m − 3 .

Fig. 2.
Environmental variables monitored at the sampling locations (shown as means of 4 m and 40 m from ditch) during GHG flux measurements averaged for each sampling occasion during year 2018 and 2019. The variables include (a) air (T a ) and soil temperature (T s5 at 5 cm and T s10 at 10 cm depth), (b) ambient photosynthetically active radiation (PAR), (c) soil moisture at 5 cm depth (SM), (d) water table level. For (c) and (d), data were grouped by ditch treatment uncleaned (UC) and cleaned (DC). Means of each sampling occasion ± standard error (SE) for both UC and DC areas (n = 10), respectively. As manual SM measurements were not conducted in 2018, SM values from the automated sensors, calibrated using 2019 data, are presented to illustrate the mean SM level near the center of the site (See Fig. 1).

and Table A.2).
Although there were about 30 % of measurement frames without significant vegetation (areal coverage <10 %) in both UC and DC areas, half of the measurement frames in UC area was comprised of substantial amounts of vegetation (areal coverage >60 %) leading to a larger range of vegetation areal coverage and greenness index particularly within the UC area (Fig. 3). In July 2019, ground vegetation became more abundant spreading across the entire site. At that time, similar in-frame mean vegetation areal coverage values were recorded for the UC area (51 ± 12 %) and DC area (46 ± 8 %), and also for the mean in-frame greenness index, UC (0.38 ± 0.01) and DC (0.38 ± 0.01) (Fig. 3 and Table A.2). Comparing in-frame vegetation with the surrounding area, we observed similar greenness indices in the surrounding UC and DC areas (UC: 0.38 ± 0.01; DC: 0.38 ± 0.01) in July 2019. Similarly, vegetation coverage was not significantly different (p > 0.05) between the surrounding UC and DC areas (UC: 54 ± 6 %; DC: 59 ± 4 %), which is consistent with in-frame vegetation coverage. There was a large spatial variation of vegetation development within the same treatment area, but we found no significant difference (p > 0.05) in vegetation coverage between sampling locations at 4 m and 40 m from the nearest ditch, for both DC treatments in both years (Table A.2).

Greenhouse gas fluxes in clear-cut area
The seasonal variations of daytime CO 2 flux partitions (NEE, GPP and R eco ) were relatively larger in 2019 than in 2018 for both UC and DC areas (Fig. 4a-c). In 2018, the average daytime NEE resulted in emissions of up to 214 and 150 mg m − 2 h − 1 of CO 2 -C at UC and DC areas, respectively. In 2019, average daytime NEE switched to net uptake with maximum values of − 36 and − 176 mg m − 2 h − 1 of CO 2 -C at the UC and DC areas, respectively (Fig. 4a). In both UC and DC areas, peak season GPP increased by two to three times from 2018 to 2019. Similarly, R eco increased by nearly-two times in both UC and DC areas from 2018 to 2019. The peaks of GPP and R eco were recorded in late-July to mid-August for both years. There was higher pronounced seasonal variability of GPP and R eco at the sampling locations with higher mean growing season fluxes (Fig. A.9). In comparison, the temporal patterns of daytime NEE showed relatively limited seasonal variation in both years. It is further noteworthy that the enhanced spatial variability, in response to vigorous but patchy ground vegetation establishment, considerably increased the uncertainty range of the mean fluxes in 2019.
During the first two years following DC, a net uptake of CH 4 of similar magnitude ranging from − 32 to − 126 µg C m − 2 h − 1 was observed across both UC and DC areas at the respective sampling occasion (Fig. 4d). Peak uptake was observed in July and August for both years, with similar ranges of spatial variation within the same treatment area for both years. Compared to CO 2 , the spatial variability was similar in both years and caused by the range of consistently higher or lower fluxes occurring throughout the year at specific chamber measurements locations (Fig. A.10). N 2 O fluxes in 2019 occurred mostly within the 10-90 percentile range of − 14 to 34 µg N m − 2 h − 1 across both UC and DC areas, respectively (Fig. 4e). In contrast to CH 4 , the spatial variability in N 2 O fluxes was caused by sporadic high and low fluxes at a given measurement location.

Biotic and abiotic controls of temporal and spatial variations of GHG fluxes in clear-cut area
The PCA for the CO 2 and CH 4 fluxes based on 2018 and 2019 data revealed two significant principal components (PCs) from each year, explaining 63.4 % and 56.6 % of the total variance for the two years, respectively (Fig. 5). In 2018, the CO 2 component fluxes (GPP, R eco and NEE) had the largest association with vegetation greenness index (g cc ) relative to other environmental variables of WTL and soil temperature. The component loadings of ditch cleaning (DC) of PC1 and PC2 aligned at the same direction with these variables indicating its negative relationships with g cc and R eco , along with positive relationships with GPP and NEE. The other variables, namely soil temperature, WTL, distance to ditch and CH 4 flux, were instead aligned in the relatively perpendicular direction at PC1 and PC2 (Fig. 5). Particularly, high loading of CH 4 flux (i.e., small uptake of CH 4 ) was associated with high WTL, low soil temperature, and greater distance to ditches. In 2019, the strong association between g cc and CO 2 component fluxes (GPP, R eco and NEE) persisted which was mostly presented in PC1, but DC became relatively independent from these variables with its negligible contribution to PC 1. Instead, DC became more associated to lower CH 4 flux (i.e., high uptake of CH 4 ) together with lower WTL and shorter distance to ditches, all of which are positioned along the similar degree of directions.
The PCA for the N 2 O fluxes also revealed two significant PCs (Fig. A.11), explaining altogether 49.6 % of the total variance (Fig. A.11). Small component loadings N 2 O fluxes in both PC1 and PC2 indicate that the other studied variables did not have strong covariations with N 2 O. N 2 O contributed slightly to PC1 that its high score represents mainly high WTL and low soil temperatures, but with minor representation to DC and distance effects.

DC treatment effects on the GHG fluxes in clear-cut area
The mixed effect models suggested that DC significantly reduced R eco in both years (p < 0.02), as well as the in-frame vegetation area cover and GPP in 2018 (p < 0.05) ( Table 2). However, NEE was not significantly affected by DC in both years (p > 0.05). The mixed effect models further suggested that CH 4 uptake at 4 m distance from both UC and DC ditches was significantly (p < 0.01) higher than at 40 m distance (Table 2). Results from mixed effect model showed no statistically significant variation in N 2 O fluxes neither to temporal variation nor to DC treatment and distance to ditch variables, while suggesting a significant effect from the treatment-distance interaction (p = 0.03), in response to the high fluxes recorded at 40 m from uncleaned ditches (Table 2).

Ditch CO 2 and CH 4 fluxes
During both 2018 and 2019, the ditch water level in May and early June was around 4 to 27 cm higher (p < 0.05) in the uncleaned than in the cleaned ditches, and all ditches dried out during June to August (Fig. 6a). The 10-90 percentile range of CO 2 emission from DC ditches was 0.2 to 150 mg C m − 2 h − 1 and higher compared to 0.1 and 110 mg C m − 2 h − 1 observed at the UC ditches during the two study years (Fig. 6b). The annual mean (median) CO 2 emission rate in 2018 was 61 ± 18 (44) and 55 ± 11 (24) mg C m − 2 h − 1 at UC and DC ditches, and in 2019 was 52 ± 12 (31) and 59 ± 14 (30) mg C m − 2 h − 1 at UC and DC ditches, respectively.
Ditch CH 4 emissions were highly variable in the UC ditches, particularly during the occurrence of higher ditch water levels during the wet spring season, with the 10-90 percentile range between − 2.9 to 109 µg C m − 2 h − 1 during the two study years (Fig. 6c). In comparison, the 10-90 percentile range of the CH 4 flux from the DC ditches ranged from − 11 to 22 µg C m − 2 h − 1 . Annual mean (median) ditch CH 4 emission rates of 168 ± 128 (0.0) and 67 ± 50 (0.6) µg C m − 2 h − 1 in 2018 and 2019 were estimated for UC ditches, in comparison to 12 ± 6 (3.5) and 2 ± 2 (0.0) µg C m − 2 h − 1 in 2018 and 2019 for DC ditches.

Total annual greenhouse gas balance
Our model estimates suggested that DC reduced the annual GHG emissions by a mean of 44 % in the second year after ditch cleaning (i.e., 2019) relative to the UC area (Table 3). In absolute terms, CO 2 contributed the most to the total GHG budget with 99 and 98 % in UC and DC clear-cut areas, respectively. CH 4 contributed 0.6 and 1.0 %, and N 2 O contributed 0.5 and 1.3 %, to the total GHG budget in UC and DC areas, respectively. Annual GPP was 25 % lower in the DC area than in the UC area. Annual R eco was 32 % lower in the DC area than in the UC area, which resulted primarily from the 34 % lower R h in the DC area. Also, R a was 25 % lower in the DC area, estimated using a range of carbon use efficiencies of 0.5. The estimate remains at 25 to 26 % under various carbon use efficiencies from 0.4 to 0.6 (Table A.3). The modelled annual CH 4 uptake in CO 2 eq ha − 1 year − 1 was slightly but significantly higher in the DC area than in the UC area. The annual emission of N 2 O in the DC area was statistically not different from that in the UC area. Loadings, representing the measured variables, are represented by filled symbols. Scores, representing the observation units, i.e., the frame locations and timing of measurements, are represented by dots (UC area) and crosses (DC area) in grey (4 m from ditch) and yellow (40 m from ditch) colour. Note that ecosystem uptake in NEE and GPP is given negative sign resulting in negative correlation to PC1 despite positive causal relation. Abbreviations represent ground vegetation areal coverage (VC), soil 5 or 10 cm-depth temperature (T s5 , T s10 ), water table level (WTL) for environmental variables (red symbols); net ecosystem exchange (NEE), ecosystem respiration (R eco ), gross primary productivity (GPP) and methane (CH 4 ) for flux variables (black symbols); and DC treatment and distance to ditch variable (blue symbols). Note that the DC variable was fitted as a binary variable for UC (0) and DC (1) treatments.

Table 2
Annual mean ± standard error (SE) and mixed effect model results for treatment effects (UC = uncleaned; DC = ditch cleaning) and distance combinations (UC4, DC4: 4 m from a ditch, and UC40, DC40: 40 m) on GHG fluxes based on observations in 2018 and 2019. Fixed factors of mixed effect models include ditch cleaning treatment (T), distance to ditch (D) and their interaction (TD). Column N refers to the sample size of each model. Note: Degrees of freedom (df 1, df 2) = 1, N-4 for mixed effect model statistics of all variables. At each sampling occasion, treatment and distance category, the number of replicates were 5. Asterisk indicates statistical significance at α = 0.05.
The annual ditch CO 2 emission was similar for UC and DC ditches. In comparison, the annual ditch CH 4 emission remained negligible for both UC and DC ditches in terms of the GHG balance. Although the absolute values for the proportion of ditch to clear-cut area fluxes were 50.2 % and 6.5 % for CO 2 and CH 4 flux, respectively, the total contribution of ditch fluxes to the area-weighed GHG balance was minor, since the areal coverage of ditches was small (2.4 %) relative to the clear-cut area.

Initial effects on spatio-temporal variations of GHG fluxes following DC
Our measurements demonstrate that among the three studied GHGs, the CO 2 flux was the dominant component, with pronounced differences between UC and DC areas observed for both GPP and R eco . The similar daytime NEE observed in UC and DC areas during both 2018 and 2019 could be because the decreased magnitude of daytime GPP at the DC area was offset by the concurrent decrease in R eco . Lower daytime GPP in the DC area in 2018 was likely associated with the delayed development of in-frame vegetation. Unfortunately, vegetation data from the surrounding area was lacking in 2018 to confirm that this pattern was a general feature and not specific to our measurement frames. However, the good agreement between in-frame and surrounding vegetation growth in 2019 suggests that in-frame vegetation also represented the general vegetation development at the site in 2018.
The mechanism of the delayed herbaceous vegetation development following DC in our study remains uncertain, however, it is to be noted that the upper part of the soil profile at the site was relatively dry with a mean WTL of − 54 to − 56 cm in all our sampling locations in the two sampling years (Table A.2) which is in contrast to previous studies in boreal mineral soil forests (e.g. Tong et al., 2022) and drained peatland forests with shallower mean WTLs of − 30 to − 40 cm (e.g. Leppä et al., 2020;Lohila et al., 2011;Ojanen and Minkkinen, 2019). One possible explanation is that the meteorological drought stress combined with the additional soil water reduction following DC might have made it more difficult for herbaceous plants species to establish in the first year. Similarly, the lower R eco in DC area in both years (2018 and 2019) was likely due to lower R a from the decreased amount of in-frame vegetation biomass, possibly in combination with subsequently decreased plant litter input in both 2018 and 2019, limiting microbial decomposition (i. e., R h ). Alternatively, the observed differences in herbaceous plant establishment might also result from other unknown, partly stochastic factors regulating seed dispersal and germination. While a thorough examination on the vegetation responses to DC was beyond the current scope, further studies are encouraged to explore in more detail how DC affects the establishment of herbaceous plants species. It is further noteworthy that DC lowered the mean WTL by <10 cm in 2019 and the DC effect on WTL in 2018 was not significant (Table A.2). Possibly, enhanced transpiration from the earlier and more developed herbaceous vegetation in the UC area might have counterbalanced the difference occurring due to the drainage effect from DC. Moreover, although the effect of DC on WTL was relatively small and occurring at a lower soil depth, the change in WTL was significantly correlated (p < 0.05) with the reduction of soil moisture in the surface layer. This indicates the potential importance of herbaceous vegetation on the water regulation along the soil profile (Ruseckas et al., 2015). Such WTL-SM-g cc interactions were also revealed by PCA results from both years which highlight the potential consequences of the vegetation and soil water dynamics on the CO 2 exchange. Hence, it is important to consider the pre-drainage WTL condition when evaluating the potential DC effect on the vegetation growth and thus CO 2 exchange dynamics.
The drainage function commonly diminishes with distance from the ditch (Prévost et al., 1997) and consequently spatial variations in DC effects on GHG fluxes may occur. Our mixed effect model results demonstrated a significant treatment effect from ditch distance on WTL at 4 m relative to 40 m from ditches, however, there was no effect of distance to ditch on CO 2 fluxes in neither DC nor UC areas within the first two years. The lack of distance to ditch effect on CO 2 fluxes is likely explained by the fact that biotic and environmental factors controlling CO 2 production and respiration, e.g., vegetation growth and soil temperature remained similar at 4 m and 40 m from ditches in our study (Table A.2). Roy et al. (2000) found that the effect of distance to ditch on ground vegetation growth was pronounced only from the third year after clear-cut, which indicates that our initial two years of measurements might not have captured the effect of distance to ditch on vegetation growth that may become evident only several years after DC.
Compared to the CO 2 fluxes, the CH 4 uptake in DC area did not deviate substantially from the UC area and the increase was only apparent in the second year (2019) following DC. This might be explained by the delayed effect of DC on lowering WTL, the latter being a major control of CH 4 fluxes (Maljanen et al., 2010;Ojanen et al., 2010) as further evident from the PCA results and additional single-factor correlation analysis (Fig. A.12), which was significant only in 2019 (p < 0.01). The net uptake of CH 4 decreased (p < 0.01) with distance from both UC and DC ditches in both growing seasons. In the DC area, the lower WTL and soil moisture observed at 4 m compared to 40 m might explain the greater CH 4 uptake near the ditch (Table 2 and  Table A.2).
Given the nutrient rich conditions (C:N ratio = 16 ± 1.3), the small N 2 O emissions were unexpected since earlier studies have shown a high potential for N 2 O emissions for drained organic soils with C:N ratios below 25 (Gundersen et al., 1998;Klemedtsson et al., 2005). The small N 2 O emissions at this site were possibly also the result of the relatively low soil moisture content associated with a low WTL compared to other studies in drained peatland forest clear-cuts (e.g. Saari et al., 2010), since suppressed denitrification was previously observed at WTL below − 30 cm (Hefting et al., 2004). Also, since the production and emission of N 2 O is commonly highly variable in both space and time as it depends on a complex series of different processes and pathways (Robertson and Tiedje, 1987;Webster and Hopkins, 1996), it is also possible that our biweekly sampling campaigns failed to capture occasional emission events (Smith and Dobbie, 2001). Thus, high spatial and temporal resolution measurements (i.e., using automated chamber or eddy covariance flux systems) are encouraged to study the response of N 2 O emissions to post-harvest DC on drained peat soils in more detail (Pihlatie et al., 2005).
Ditch emissions of CO 2 and CH 4 : were higher and lower, respectively, compared to the range (CO 2 : 21-46 mg C m − 2 h − 1 ; CH 4 188-6838 µg C m − 2 h − 1 ) reported previously from ditches in other drained peatland forests (Ball et al., 2007;Peacock et al., 2021c;von Arnold et al., 2005). In comparison to previously studied ditches in drained peatland forest sites (e.g. Minkkinen et al., 1997;Peacock et al., 2017), the conditions in both UC and DC ditches were relatively dry which may have enhanced soil aeration, resulting in increased CO 2 respiration and aerobic CH 4 oxidation while suppressing anaerobic CH 4 production (Nykänen et al. 1998). It is however noteworthy that the CH 4 emission spikes occurring in the UC ditches during wet spring periods were not observed in the DC ditches. The CH 4 emission spikes in the UC ditches might have been facilitated by the water-logged conditions, which create a favourable anaerobic environment for CH 4 production. Furthermore, the presence of herbaceous species in the UC ditches may have enhanced root-derived substrate supply thereby supporting greater CH 4 production (Zhang et al., 2002).

Total GHG balance following DC
Our study revealed a strong reduction of annual net GHG emissions by almost half in the second year following DC compared to the UC treatment. Given the lack of extensive pre-treatment data, we need to caution and acknowledge that it remains somewhat elusive whether this observation is of causal or correlative nature. However, the lack of differences in pre-treatment WTL and soil chemistry (or any other of the measured environmental properties) among treatment plots strongly suggests that the altered GHG balance is a result of DC. We further acknowledge that the absolute annual GHG budget sums are somewhat uncertain given the need for extrapolating fluxes during the winter months. However, relatively low baseline emissions can be expected during this period, thus resulting in only a limited contribution to the annual sums. Furthermore, given similar model parameter uncertainties and by using the same approach for extrapolating winter fluxes for both treatments, it is unlikely that the model extrapolation affected the results on the DC treatment effect on the annual GHG balances.
The main reason for the observed reduction in the total annual GHG balance following DC was primarily the lower soil CO 2 emissions, being Table 3 Model estimates for the total annual GHG balances in 2019 based on exchanges of methane (CH 4 ), nitrous oxide (N 2 O) and carbon dioxide (CO 2 ) fluxes, the latter including the net CO 2 exchange (NEE) and its component fluxes of gross primary productivity (GPP) and ecosystem respiration (R eco ) which comprises heterotrophic respiration (R h ) and autotrophic respiration (R a ) in uncleaned (UC) and cleaned (DC) ditch treatments. Values are in the unit of t CO 2 eq ha − 1 year − 1 using global warming potentials of 34 and 298 for CH 4 and N 2 O over a 100 year timeframe, respectively (IPCC, 2013). Numbers are represented with ± standard deviation (SD) based on Monte Carlo uncertainty estimates. a The area-weighted sum includes the relative contribution of fluxes from ditch (width = 1 m) and two distances (area of 0-22 m from ditch based on 4 m flux and remaining area based on 40 m flux) from ditch, in a ditch network at 185 m intervals. b Photosynthetic CO 2 uptake and plant R a from occasional plants and mosses growing in the ditches were not accounted for and assumed to be zero in the ditch net CO 2 exchange. c Carbon use efficiency (CUE) of 0.5 is assumed to derive R a = GPP × 0.5 and (Waring et al., 1998;Gifford, 2003). R a and R h estimated using CUE = 0.4 and 0.6 are shown in Table A.3. d Ditch N 2 O fluxes were not determined. e The total GHG balance is the sum of the area-weighted (combining clear-cut areas and ditches) exchanges of CO 2 (i.e. NEE), CH 4 and N 2 O, assuming a constant ditch width of 1 m at 185 m intervals within the ditch network.
the dominating component of the total GHG balance, from the DC area. The dominating contribution of the CO 2 exchange to the ecosystem GHG balance (over a 100 year time frame) is in line with earlier studies of the GHG balance in clear-cut sites on drained peatlands (e.g. Korkiakoski et al., 2019;Vestin et al., 2020). Our mass balance approach assuming a plant CUE of 0.5 at the annual scale (to derive R a from GPP, and R h from the difference of R eco and R a ) suggests that a considerable decrease in R h (by 34 %) was likely the main reason for the decrease in the total GHG balances after DC. It is noteworthy that we arrived at the same results when assuming a plant CUE of 0.4 or 0.6 (Table A.3), the latter spanning the commonly reported range in CUE (DeLucia et al., 2007). A decrease in R h following WTL draw-down seems at first counter-intuitive since numerous studies previously suggested that drainage of wet peat soils enhances microbial decomposition of soil organic matter (Maljanen et al., 2010;Ojanen et al., 2013;van Huissteden et al., 2006). However, since our study site appeared relatively dry with the lowest WTL reaching − 107 cm during the summer 2019, it is possible that the increased dryness following DC at our specific site might have led to drought-induced inhibition of microbial activity and reduced R h (Drzymulska, 2016, Manzoni et al 2012. In addition, the extensive cover of herbaceous plants species (e.g., Chelidonium majus) established within the UC area might have supported heterotrophic decomposition by providing enhanced input of easily decomposable organic matter (Thormann et al., 2001). Altogether, this indicates that the enhanced drought stress after DC activities might have suppressed rather than stimulated R h and net CO 2 emissions at our site, which was characterized by relatively low pre-DC WTL and experienced a meteorological drought in the study year 2018. These findings are in contrast to DC effects at wetter sites where a reduction in WTL following DC might enhance soil mineralization and CO 2 emission rates (Maljanen et al., 2010;Ojanen et al., 2013;van Huissteden et al., 2006). Thus, DC effects on the CO 2 balance might vary across different sites in dependence of the combined effects from the effectiveness of the pre-DC ditch drainage function, weather patterns and site hydrological conditions. Compared to the CO 2 fluxes, the contribution of soil CH 4 and N 2 O fluxes to the total GHG balance remained limited despite their considerably higher warming potential over a 100-year time frame. Even if considering the warming potential for CH 4 of 86 over only a 20-year time frame (Myhre et al., 2013), its contribution was < 2 %. The negligible contribution of CH 4 fluxes to the total annual GHG budget noted in our study is in line with previous studies on recent forest clearcuts (e.g. Korkiakoski et al., 2019;Vestin et al., 2020). However, under wetter conditions the response of CH 4 fluxes after clear-cutting might be higher (Bradford et al., 2000) compared to the dry conditions at the studied site. Furthermore, due to the slow recovery of the methanotrophs to water table fluctuations (Adamsen and King, 1993), the contribution of CH 4 uptake might become more pronounced in later years (Koschorreck and Conrad, 1993;Whalen and Reeburgh, 1996;Gulledge and Schimel, 1998). Similarly, the small contribution of N 2 O to the total GHG balance (<1.5 %) might be the result of the dry conditions at this site whereas in wetter (and nutrient-rich) sites a greater contribution and response to DC could be expected (Martikainen et al., 1993;Regina et al., 1996;Regina et al., 1998), but more spatio-temporal replication is needed to capture the influence of potential spatial and temporal emission hotspots on the total GHG balance.
Previous studies suggested that GHG emissions from ditches may considerably affect the ecosystem GHG balance, predominantly through high CH 4 emissions (Hyvönen et al., 2013;Minkkinen and Laine, 2006;Peacock et al., 2017;Peacock et al., 2021c;Schrier-Uijl et al., 2010;Teh et al., 2011). At this site, however, the lack of high GHG fluxes from ditches in combination with the small ditch area (2.4 %) resulted in minor area-weighted contributions from ditch CO 2 (1.5 %) and CH 4 (<1%) emissions to the total GHG budget. The minor contribution of ditch CH 4 emissions to the total GHG balance was likely because the ditches at our dry site frequently dried out during the study period, which resulted even in occasional CH 4 uptake.
To our knowledge, there is only one published study to date that explored DC effects on the forest GHG balance which was carried out in a forest clear-cut on wet mineral soil in boreal Sweden (Tong et al., 2022). This study reported insignificant effects of DC on C and GHG balances and vegetation development, despite changes in WTL due to DC. Given the contrasting results from this northern site compared to ours, additional studies across a range of varying site characteristics are therefore urgently needed to obtain a more generalized understanding of DC impacts on forest GHG dynamics. Furthermore, while our study only addressed the initial responses, additional effects from the enhanced tree growth might further modify the DC impact on the forest GHG balance over an entire stand rotation. These impacts may include further soil water reduction due to increased evapotranspiration, increased photosynthetic CO 2 uptake by trees but also greater canopy shading effects on soil temperature and moisture levels. In addition, the recently cleaned ditches might degrade over time and eventually lose their drainage function which will feedback on the soil water dynamics and hence GHG dynamics. Thus, more empirical data are needed to fully understand the long-term DC effects on the forest ecosystem GHG balance. Given the steady increase of DC activities within Fennoscandia, a better understanding of its effects on the forest ecosystem-atmosphere exchange of GHGs is urgently needed to support the development of appropriate forest management strategies to mitigate climate change.

Conclusions
We examined post-harvest ditch cleaning (DC) effects on CO 2 , CH 4 and N 2 O fluxes from clear-cut area and ditches in a drained fertile peatland forest site in hemiboreal Sweden over two years after clearcutting and DC. Based on our findings we conclude that: 1. The effect of distance to ditch was insignificant for CO 2 and N 2 O fluxes while being small but significant for CH 4 fluxes during these two initial study years. Particularly for the CO 2 exchange, the lack of a clear spatial response to ditch location was likely due to overshadowing effects from the high spatial variability in ground vegetation establishment. 2. Soil water dynamics (i.e., WTL and soil moisture) and ground vegetation coverage were identified as the main controls on the spatial variations of measured CH 4 uptake and CO 2 component fluxes (R eco and GPP), respectively. Specifically, lower soil water content and the delayed vegetation development in the DC area coincided with larger CH 4 uptake and smaller CO 2 component fluxes (i.e., production and respiration), respectively, relative to the UC area. 3. Model extrapolations suggest that, in the second year following clear-cutting, the drained peat soil acted as considerably GHG source, however, with 44 % lower emissions from the DC area (27.8 ± 10.3 t CO 2 -eq ha − 1 year − 1 ) compared to the UC area (49.4 ± 17.0 t CO 2 -eq ha − 1 year − 1 ). While direct evidence is lacking, we propose a decrease in R h (by 34 %) due to enhanced soil water stress at the relatively dry study site as a potential reason for the reduction of total GHG emission in the DC area, relative to the UC area. 4. Overall the ecosystem GHG balance was dominated by the CO 2 exchange. The contribution from fluxes of CH 4 and N 2 O from the clearcut area as well as CH 4 fluxes from ditches (N 2 O not sampled) were negligible.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.