2.1 The site and management

The site was located in southern Finland (60.22° N, 24.78° E, 110 m a.s.l.), and it has been in cultivation at least since the 19th century. The field has been in a crop rotation with cereals and grass during the latest decades. The climate is boreal humid, with long-term (1991–2021) annual mean temperature of 5.2 °C and precipitation of 621 mm (Jokinen et al., 2021). The sum of annual global radiation is 3358 MJ m⁻² and total sunshine duration 1699 h. Typically, the soil is frozen and has a snow cover from December to March–April. The field was a highly decomposed fen, with peat depth ranging from 0.8 to over 2 m. Organic carbon content was 29 % and pH 5.4 in the surface layer (0–20 cm) (Table 1). The original subsurface drainage system with tile drains was replaced by modern plastic pipes surrounded by gravel in the 1960s. The distance between the pipes was 18 m until 1979, when it was changed to 9 m. The drainage depth was 60–80 cm, and a control well was installed prior to the experiment to restrict water outflow and raise the groundwater table. The adjustable tube inside the well was set to a position that let the water out when the water table reached 20 cm depth below the soil surface.

Table 1Soil properties (± standard deviation) in the 0–20 cm layer in 2021.

Download Print Version | Download XLSX

The site was established in 2018, and it consists of 12 experimental plots (9×6 m) in four blocks (see the key figure). Four replicate plots with either a grass mixture for forage (sown with Poa trivialis and Festuca pratensis, replanted in 2019 and 2021 with Phleum pratense, Festuca pratensis, Lolium multiflorum and Poa pratensis), bog bilberry (Vaccinium uliginosum; a.k.a. bog blueberry or bog whortleberry) or willow variety Klara (hybrid of Salix schwerinii L. `Amgunskaja' × Salix viminalis L. `Ivar') were randomly assigned within the four blocks. The grass was seeded and bilberry seedlings and willow saplings planted in June 2018 (Table S1 in the Supplement). The bilberry did not grow roots, and those plots were left to develop to the set-aside type during the following years; thus, we named this treatment as set-aside. The number of species in all 12 plots was determined once in the summer of 2021.

2.2 Ancillary measurements

Biomass growth of willow was monitored by cutting three willow individuals from each plot for determining the above-ground biomass each June. The leaves and stem with branches were separated and weighed to determine fresh biomass. The woody biomass was cut in 10 cm pieces and dried at 65 °C for 2 weeks. The root biomass around one of the monitored plants per plot was determined by taking $\mathrm{50}\times \mathrm{80}\times \mathrm{20}$ cm peat samples from three layers: 0–20, 20–40 and 40–60 cm once per year. Visible large (>2 mm) and fine roots were manually separated from the peat, dried and weighed. For determining fine roots, the peat samples were mixed, and a 1 kg subsample was taken. Annual growth in stem, stool and coarse roots was calculated by subtracting the value from the previous year. Annual turnover rate of fine roots was assumed to be 3 times the biomass of fine roots as in Pacaldo et al. (2014). For example, biomass increment in 2019 was calculated with the following equation:

where F₁₉ is foliage in 2019, S₁₉ and S₂₀ are stems in 2019 and 2020, St₁₉ and St₂₀ are stools in 2019 and 2020, Cr₁₉ and Cr₂₀ are coarse roots in 2019 and 2020 and Fr₁₉ is fine roots in 2019. Subsamples were taken for determining the C content of the dried biomass in 2019 and 2020, and the mean values were used for the following years. The yield per hectare was estimated to be the weight of 25 000 individuals based on 80 cm ×50 cm spacing.

Soil temperature was measured first at a depth of 10 cm (but at a depth of 5 cm from May 2020 on to achieve better response of CO₂ to air temperature) in each treatment with ElcoLog sensors (Elcoplast Oy, Tampere, Finland). The sampling rate was 1 h in summer and 2.5 h in winter. The air temperature, precipitation and radiation data were taken from the Jokioinen weather station of the Finnish Meteorological Institute (FMI, 2024, CC BY 4.0) located about 10 km from the site. Continuous photosynthetically active radiation (PAR) data were produced with global radiation data from FMI and corrected using the ratio of 2.04 for global radiation and PAR (Meek et al., 1984).

WTD was measured from monitoring pipes at the corners of the site at the time of the opaque chamber measurements until 2021, when monitoring pipes were also installed in the centre of each plot. During the summers of 2021 and 2022, there were also HOBO water level data loggers (Onset, Bourne, MA, USA) in each plot for continuous water level monitoring with a sampling rate of 1 h. In winter, when the loggers were not used, WTD was measured manually from monitoring pipes when the water was not frozen.

The leaf area index (LAI) was measured at the same time as the transparent chamber measurements with a portable LAI metre (SunScan; Delta-T Devices Ltd, Cambridge, United Kingdom). LAI values of >3 were set to 3 as they were assumed to not affect photosynthesis due to the saturation of the reflectance (Aparicio et al., 2000). When harvesting the grass plots, the previous measured LAI value was extrapolated to the moment just before harvesting, after which the LAI value was set to 1 as measured. In 2022, we measured green canopy cover with the Canopeo app (Patrignani and Ochsner, 2015) instead of LAI. Based on our experiences and due to the operation and physical design of the LAI device, it did not provide as comprehensive of a picture of the biomass inside the gas measurement collar as Canopeo. The vegetation index was found to be faster to measure and less dependent on the ambient light conditions than the light interception method (Shepherd et al., 2018). LAI was indexed by dividing by the maximum value of 3 and green canopy cover by 100 (values from 0 % to 100 %) so that the generated vegetation index range was 0–1. The vegetation index was set to 0 from the end of November until mid-April, when snow and frost covered the ground or no green vegetation was present.

Soil samples for analysing the soil properties were taken first in October 2018, and another sampling was conducted in June 2021 with additional analysis. As there were no significant changes in the soil properties between these samplings, we only present the results of the second soil sampling in Table 1. The samples were taken from the 0–20 cm layer using a soil corer with a diameter of 3 cm. Approximately 20 subsamples were pooled to make a composite sample that was air-dried and sieved (2 mm) for the chemical analyses. Soil core samples for dry bulk density and porosity (diameter of 5 cm) were taken from the surface layer (0–17.5 cm) of each plot in October 2020 using the Kopec corer, and the samples were dried at 37 °C for a week. Soil acidity was determined using the ISO 10390 method. Nutrient content was analysed as described in Vuorinen and Mäkitie (1955). Soil carbon and nitrogen were determined using the dry combustion method (Leco TruMac CN; LECO Corporation, MI, USA).

2.3 GHG measurements

Dark respiration of the plants together with soil respiration (ecosystem respiration) and fluxes of N₂O and CH₄ were measured using opaque chambers biweekly or once per month in the winter between March 2019 and March 2023. In each plot, a 60 cm×60 cm steel collar was installed at a depth of 10–15 cm. The location of the collars was 1 m from the short edge of the plot and 3 m from the edges of the adjacent plots. An aluminium chamber (height of 40 cm) mounted at the top of the collar was sealed with water in the groove of the upper edge of the collar. In the winter, NaCl was added to the water to avoid ice formation. The clear aluminium surface effectively reflected light and kept the temperature change moderate inside the chamber. The measurements were done during the daytime between 10:00 EET and 14:00 EET approximately every 2 weeks in summertime and monthly in the winter. The chambers were closed for 30 min, and four 20 mL gas samples were taken with a 60 mL plastic syringe to pre-evacuated vials (Exetainer, Labco Limited, UK) at 10 min intervals starting immediately after closing. Prior to sampling, the syringe was pumped five times to mix the air in the chamber. The samples were analysed with a gas chromatograph (Agilent 7890; Agilent Technologies, Inc., Wilmington, DE, USA), equipped with a flame ioniser and electron capture detectors and a nickel catalyst for converting CO₂ to CH₄. The gas chromatograph had a 2 mL sample loop and a backflush system for separating water from the sample and flushing the precolumn between the runs. The precolumn and analytical columns consisted of 1.8 and 3 m long steel columns, respectively, and were packed with an 80/100 mesh HayeSep Q (Supelco Inc., Bellefonte, PA, USA). Nitrogen was used as the carrier gas, and a standard gas mixture of known concentration of CO₂, N₂O and CH₄ was used for a calibration curve with seven concentration points. An autosampler (222 XL Liquid Handler; Gilson Medical Electronics, France) fed the samples to the loop of the gas chromatograph.

Net ecosystem exchange (NEE) including photosynthesis and respiration of the soil and plants was measured approximately every 2 weeks during the growing season using a transparent chamber ($\mathrm{60}\times \mathrm{60}\times \mathrm{60}$ cm) made of polycarbonate plexiglass (1 mm, light transmission of 95 %). The chamber was equipped with a Vaisala GMP-343 probe for CO₂ measurement and a temperature and humidity sensor (Vaisala Oy, Vantaa, Finland) and two fans for mixing the air during the measurement. PAR was measured with a LI-190 quantum PAR sensor (LI-COR, Lincoln, NE, USA) inside the chamber. Four measurements with different amounts of light entering were taken from each plot on each measurement day in order to cover a large range of light conditions and facilitate the gap filling by modelling. One or two layers of a white fabric shroud and one blackout curtain were used to acquire measurement results in different light conditions (approximately 100 %, 50 %, 25 % and 0 % of ambient radiation). The measurement with 0 % radiation gave the estimate for ecosystem respiration (ER). The measurements were done in the same collars as the opaque chamber measurements. Each measurement took 1 min with a 5 s sampling rate or 2 min in early or late growing season when the change in flux was minor. The chamber was flushed after each measurement to reconstitute ambient CO₂ and air humidity contents. After closing the chamber, a lag time of 10 s was applied to exclude the time when the flux was not yet stabilised. Clear-sky conditions were preferred to avoid problems related to changing cloud cover and to achieve the widest possible range of available light. The temperature change inside the chamber was less than 1.5°, which was also used as a criterion for data filtering.

The change in CO₂ concentration during the chamber enclosure was assumed to be linear. The measurement results of CO₂ in parts per million (ppm) were converted to g m⁻² h⁻¹ by the ideal gas law using measured temperature inside the chamber. If the flux was not yet stabilised at the beginning (first four data points) of the measurement, outliers were defined with the MATLAB isoutlier command, resulting in the removal of 210 of total 23 066 data points in 1564 flux measurements.

If the snow cover was thicker than 20 cm, a concentration snow gradient method as in Maljanen et al. (2003) was used to determine the GHG fluxes. A probe made of a steel pipe (Ø 3 mm), with a three-way valve and a plastic syringe, was used to sample 15 mL of air just above the snow cover, at the bottom of the snow cover and at every 10 cm in between in three replicate locations per plot. The gas was stored in the pre-evacuated vials, and the concentrations were determined by gas chromatography.

Measurements for bare soil respiration were made in unvegetated subplots in July 2019–December 2022. For willow, the large 60×60 cm frames were used, but for forage and set-aside, we installed one sheet-metal air ventilation pipe 27 cm in diameter and 30 cm in length to the depth of 5–10 cm next to the opaque chamber collars in the eight plots of grass and set-aside. All green vegetation within the chamber area was removed, and root growth was limited by cutting around the chamber occasionally with a knife. For the measurements, the cylinders were closed, with a cover equipped with a CO₂ sensor (GMP-343; Vaisala Oy, Vantaa, Finland) and a small fan. One measurement lasted for 1 min with a 5 s sampling rate. Measurements were taken about once every week or two and more frequently in summer than in winter. In winters of 2021–2022 and 2022–2023, this method was not used due to too high a snow depth, but measurements with the snow gradient method were utilised (Maljanen et al., 2003).

2.4 Flux modelling

Gross photosynthesis (GP) can be determined as the difference between NEE and ER. Instantaneous GP was estimated for each measurement occasion by Eq. (2):

where the full darkened transparent chamber measurement result (ER) is subtracted from the light-dependent flux (NEE) measured during the same day. Thus, we follow the sign convention with positive ER and negative GP values.

The gaps in GP and ER data between the measurement occasions were predicted using hourly time series of the ancillary data. Hourly time points for vegetation index, WTD and soil temperature were acquired from the measured values by linear interpolation. Gaps in soil temperature were filled with the modified soil temperature model (Zheng et al., 1993) using the air temperature. Air temperature and PAR were assumed to be the same for all plots, whereas we used a plot-specific vegetation index and the soil temperature from the certain treatment. Hourly ER and GP were modelled using nonlinear regression (fitnlm function in MATLAB) for all eight plots in forage and set-aside treatments. Empirical models were used for ER as in Lohila et al. (2003) and for GP as in Kandel et al. (2013). Instead of the phytomass indices used in the above publications, we used the vegetation index formed according to the LAI and Canopeo measurements (index from 0 to 1) to describe the stage of the crop growth.

We used the following equation first defined by Long and Hällgren (1993) for GP to estimate empirical coefficients (A_max and k):

where PAR is the measured photosynthetically active radiation (µmol m⁻² s⁻¹), VI is the vegetation index, A_max is the asymptotic maximum (g CO₂ m⁻² h⁻¹), and k is a half-saturation value (µmol m⁻² s⁻¹). T_Scale represents the temperature sensitivity of photosynthesis and follows the equation presented by Raich et al. (1991):

where T is the measured temperature, T_min is the photosynthetically active minimum temperature of −2 °C, T_max is the maximum of 40 °C and the optimum is 20 °C as in Kandel et al. (2013).

ER was estimated using data from the opaque and fully darkened transparent chambers. The empirical coefficients (R0_s, R0_p, E0_s and b) were estimated with a nonlinear regression model similarly to the case of GP. Annual fluxes were computed as sum of the hourly fluxes with a trapezoidal method (trapz function in MATLAB v2019b). ER consists of autotrophic respiration (R_auto), i.e. plant respiration, and heterotrophic respiration (R_hetero), i.e. soil respiration (Lloyd and Taylor, 1994), with the extension of WTD as in Karki et al. (2014):

where T_soil is the measured soil temperature, VI is the vegetation index, T_air is the measured air temperature, R0_s is soil respiration at the reference temperature of 10 °C (g CO₂ m⁻² h⁻¹), R0_p is plant respiration at the reference temperature at 10 °C (g CO₂ m⁻² h⁻¹), b is the effect of WTD, E0_s is ecosystem sensitivity and b_d is the temperature dependence of dark respiration set to 5000 as in Lohila et al. (2003). Bare soil respiration was estimated like ER but using only Eq. (6). The estimated parameters R0s, Es and b and model R² are shown in Table S3.

2.5 Data processing and analysis

For the transparent chamber measurements, the criterion R²>0.9 for the fitted linear assumption of flux measurements would exclude a large amount of data, especially with a small change in CO₂, leading to a biased dataset. Therefore, we decided to add the criterion S_xy<2.3 g CO₂ m⁻² h⁻¹ for the dataset like in Kutzbach et al. (2007) (S_xy is the standard deviation of the residuals, and 2.3 g m⁻² h⁻¹ is the 95th percentile of measurements). This procedure resulted in the removal of 59 values out of total 1467 measurements. In the modelling phase, fitted values were examined, and outliers were removed to avoid distortion. Outliers were defined as observations with an absolute value of standardised residuals greater than 3. In 2019, 3 out of 260 GP values and 3 out of 243 ER values were removed. In 2020, none of the 200 GP values and 2 of 230 ER values were removed. Overall, 2 out of 365 GP values and 4 out of 247 ER values were removed in 2021, and 12 out of 583 GP values and 2 out of 323 ER values were removed in 2022. The model's estimated parameters A_max, k of GP, R0_s, R0_p, E_s and b of ER and model correlations are shown in Table S2. The measured versus model-predicted values of GP and ER are shown by treatments and years in Fig. S1 in the Supplement.

For bare soil respiration measurements in set-aside and forage, the same criteria were used as for the transparent chamber (R²>0.9 and S_xy<95 %), leading to a removal of 12 values of the total 601. In bare soil measurements in mid-summer 2022, 24 flux measurements occurred (nine values in one plot and none–five in others) of the total 147 values, which were inexplicably high (3–18 g CO₂ m⁻² h⁻¹). Values were of the same magnitude as values measured immediately after ploughing in Honkanen et al. (2023). We decided to remove these values as outliers to avoid model distortion. Soil respiration of willow was defined with the opaque chamber method, and such outliers did not occur in these measurements. In the modelling phase, outliers defined as observations with an absolute value of standardised residuals greater than 3 were removed, resulting in the removal of 13 measurements of the total 984 measurements (including all plots).

A linear regression model was fitted to calculate gas concentrations, and the ideal gas law was used to solve the flux rate for every enclosure of the opaque chambers. Nonlinear responses of CO₂ indicated a leaking chamber or other problem in the measurement, and thus, if the R² of CO₂ was less than 0.9, the results of CH₄ and N₂O were also discarded. In addition, sudden variations in CH₄ fluxes due to ebullition were filtered by selecting only flux rates with an intercept between 1.5 and 2.4 ppm. These criteria resulted to 176, 117 and 118 discarded values out of 1044 in the case of CH₄, CO₂ and N₂O, respectively. All data cleaning and processing was done with MATLAB (v2019b; The Math Works, Inc.).

2.6 GHG balance

The annual net ecosystem carbon balance was constructed as the sum of the hourly values of NEE and yield data for each year in the case of forage and set-aside treatments. Modelling was used to fill the gaps between the measurement occasions to create a continuous series of hourly values. For willow, annual estimates of carbon loss in soil respiration were available from the chamber measurements from the unvegetated frames, but as the willows were too high for the chambers, their net production had to be estimated based on biomass accumulation during 4 years (Pacaldo et al., 2014). The presented net ecosystem carbon balance of willow is thus the sum of average annual CO₂-C from soil respiration and average annual amount of carbon bound in the biomass during the first 4 cultivation years. The cumulative annual fluxes of CH₄ and N₂O for each management practice were calculated by interpolating the emissions between consecutive sampling days. Global warming potentials 27 and 273 were used for CH₄ and N₂O, respectively, to convert the results to CO₂ equivalent for the total GHG balance (Forster et al., 2021).

2.7 Statistical analyses

Linear mixed models were used to find variables explaining variation in the gas fluxes. Crop, year, WTD and all their interactions were denoted as fixed effects. Block and block–year interaction were assumed to be independent and normally distributed random effects. The most suitable covariance structure was chosen using Akaike's information criterion (AIC). The models were fitted using the residual maximum likelihood (REML) method, and degrees of freedom were estimated using the Kenward–Roger method. The residuals were plotted against the fitted values, and the normality of the residuals was checked using box plots. The data were log-transformed when needed to normalise the distributions. The method of Tukey–Kramer was used for all pairwise comparisons of means with a significance level of 0.05. After the first model run with all relevant variables, the non-significant variables were removed one by one to find the most relevant effects. All statistical analyses were performed using the SAS Enterprise Guide v7.1 (SAS Institute Inc., Cary, NC, USA).

Figure 1Daily mean of soil temperature and precipitation, photosynthetically active radiation (PAR), vegetation index, water table depth (site mean ± SD), gross photosynthesis, measured (dots) and model-predicted (line) ecosystem respiration (soil respiration for willow), and net ecosystem exchange. Annual modelling periods (April–March) are marked by a light-grey or white background.

Download

3.1 Climate and site variables

Annual mean temperature was 6.9, 6.0, 5.8 and 5.8 °C and annual precipitation 750, 600, 660 and 546 mm in 2019–2022, respectively. The number of days with snow cover on the soil within each modelling year (April to March) was 13, 81, 108 and 118, respectively. The annual mean temperature during the study years was higher than the long-term average of 5.2 °C in 1991–2020 (Jokinen et al., 2021). Two study years exhibited lower and two higher annual precipitation as compared to the long-term mean of 621 mm. The WTD showed an increasing trend in time and high within-year variation (Fig. 1). The average WTD was −54, −41, −39 and −27 cm in 2019–2022, respectively. WTD varied from −89 to −4, −77 to 2, −120 to 1.4 and −100 to 1.8 cm in 2019–2022, respectively.

The forage yields were 6.3±0.9, 8.9±0.7, 11±0.8 and 9.4±0.9 Mg dry matter ha⁻¹ in 2019–2022, respectively. There were two harvests in 2020 and three in the other years. The plots were dominated by Phleum pratense and Festuca pratensis in 2021. The dry mass yields of willow were 30±14 and 73±28 Mg dry matter ha⁻¹ in the harvests of February 2021 and 2023. Most of the C accumulation occurred in the stem (59 %) followed by stool (25 %), roots (9 %) and foliage (7 %) (Table 2). Vegetation of the set-aside plots in 2021 was dominated by wild plants belonging to Families Asteraceae, Cichoriaceae and Caryophyllaceae. Bog bilberry covered 1 % or less on each of the four replicate plots. The set-aside vegetation had the highest species diversity involving 19 vascular plants compared to 12 at willow plots and 9 at forage plots, with the latter two including crop plants.

Table 2Four-year cumulative carbon balance of willow (± standard deviation). Negative sign indicates sequestered carbon and positive sign released carbon to the atmosphere.

Download Print Version | Download XLSX

3.2 Carbon balance

Model-predicted maximum hourly GP was −0.7, −3.2, −4.3 and −4.8 g CO₂ m⁻² h⁻¹ in the set-aside plots in 2019–2020, and −3.9, −5.9, −4.3 and −4.8 g CO₂ m⁻² h⁻¹ in the forage plots in 2019–2022, respectively (Fig. S1). The maximum measured GP value was −1.1, −2.4, −3.4 and −4.5 g CO₂ m⁻² h⁻¹ for set-aside and −3.4, −6.2, −4.8 and −4.2 g CO₂ m⁻² h⁻¹ for forage in 2019–2020, respectively. Annual values of GP varied from −9.3 to −12 Mg CO₂-C ha⁻¹ yr⁻¹ in the forage and from −1.5 to −10 Mg CO₂-C ha⁻¹ yr⁻¹ in the set-aside treatment (Table 3). The variables initially included in the analysis were annual mean WTD, crop type and year as the main effects and all their interactions. Finally, only crop, year and their interaction were left in the model (Table S4). Inclusion of WTD did not result in a meaningful estimate as both the productivity of especially the set-aside established in 2019 and the WTD increased with years, and the observed change in GP did not, in reality, increase with WTD but with time. Forage and set-aside treatments differed significantly (p<0.001), there was an increasing trend in productivity from 2019–2022 (p=0.0009), and the differences between the crop types were the highest in 2019–2021 (p<0.001).

Table 3The estimated annual sums (± standard deviation) of gross photosynthesis (GP), ecosystem respiration (ER), net ecosystem exchange (NEE), carbon exported in the harvested yield, net ecosystem carbon balance (NECB), N₂O and CH₄ effluxes and the total emissions (global warming potential of 100 years; GWP-100), with either NEE or NECB representing CO₂ emissions in the forage and set-aside plots, and selected data for willow. Significant differences (p<0.05) between treatments within a year are denoted by different letters (n=4).

^a All components of the carbon balance are not available for willow; see Sect. 2.6. ^b NEE representing CO₂. ^c NEBC representing CO₂.

Download Print Version | Download XLSX

Modelled maximum hourly ER was 2.4, 2.3, 3.0 and 4.7 g CO₂ m⁻² h⁻¹ in the set-aside and 3.5, 3.4, 4.8 and 4.5 g CO₂-C m⁻² h⁻¹ in the forage plots (Fig. S1). Measured maximum ER with the opaque chamber method was 1.7, 2.0, 2.9 and 4.6 g CO₂-C m⁻² h⁻¹ for the set-aside and 3.5, 4.2, 4.0 and 4.8 g CO₂-C m⁻² h⁻¹ for forage plots. Annual ER varied from 14 to 19 Mg CO₂-C ha⁻¹ yr⁻¹ in the forage and from 8 to 17 Mg CO₂-C ha⁻¹ yr⁻¹ in the set-aside treatment (Table 3). Variables initially included in the analysis were annual mean WTD, crop type and year as the main effects as well as their interactions. WTD did not explain the variation in the annual ER estimate well for likely the same reason as for GP since plant respiration is related to the biomass of the vegetation, which increased during the experimental years. The best model was based on the crop type and year as the main effects and their interactions (Table S4). In this model, crop type explained ER well, ER increased in time and ER increased between years more in forage than set-aside.

Figure 2Plot-wise mean annual fluxes of CH₄ (a), N₂O (b) and soil respiration (c) (CO₂ eq.) as related to the mean annual WTD.

Download

Hourly model-predicted NEE varied from −2.9 to 3.7 g CO₂ m⁻² h⁻¹ in the set-aside and from −4.6 to 3.7 g CO₂ m⁻² h⁻¹ in the forage treatment (results not shown). There were 29, 28, 16 and 47 days annually with negative daily NEE in the forage plots during the study years, respectively, and fewer such days (1, 0, 9 and 6) in the set-aside plots in 2019–2022 (Fig. 1). The cumulative annual balance ranged from 5.1 to 8.0 Mg CO₂-C ha⁻¹ yr⁻¹ in the forage and from 6.8 to 7.1 Mg CO₂-C ha⁻¹ yr⁻¹ in the set-aside treatment (Table 3), and the treatments did not statistically differ. The net ecosystem carbon balance (NECB) that accounts the amount of carbon exported in the harvested yield varied from 5.5 to 13.5 Mg CO₂-C ha⁻¹ yr⁻¹ in the forage treatment and was equal to NEE in the set-aside treatment (Table 3). The NECB values differed statistically between the forage and set-aside treatments across all years (p>0.001), and inclusion of additional effects in the analysis did not improve the model (Table S4).

Annual sum of respiration varied from 8 to 15 Mg CO₂-C ha⁻¹ yr⁻¹ in the different treatments and years (Table 3). The proportion of soil respiration of the total ecosystem respiration varied from 45 to 90 % in the forage plots and from 85 % to 100 % in the set-aside plots in 2019–2022. In the set-aside plots, estimated annual bare soil respiration exceeded the estimated ER in all plots in 2019, two plots in 2020 and one plot in 2021 and 2022, and those values were not used in the above calculation, and thus it is assumed that total respiration constituted of only soil respiration in 2019. Annual cumulative soil respiration was explained by WTD (Fig. 2; p=0.053) and crop type (p=0.033), so forage and set-aside treatments were significantly different in the whole dataset and respiration increased in the order forage < willow < set-aside (Fig. S4; Table S4). Plots of the bare soil respiration in relation to WTD and temperature show that there is a clear trend of decreasing respiration with rising WTD (Fig. S2). Three individual curves indicate a contrasting trend, but these three estimations are based on a small number of measurement results. Based on all annual estimates of soil respiration, a 0.1 m rise in WTD reduces respiration by 0.87 Mg CO₂-C ha⁻¹ yr⁻¹.

The cumulative total amount of C in the above- and below-ground willow biomass was 86.5 Mg C ha⁻¹ during the 4 study years (Table 2). About 40 % of the carbon in the biomass was left at the site after harvest, and soil respiration amounted to 43.5 Mg ha⁻¹, leading to a strongly negative cumulative NEE of −43 Mg ha⁻¹. Carbon export in the harvest changed the net balance to a net loss of 7.6 Mg, corresponding to an average annual CO₂ rate of 7 Mg of CO₂.

Figure 3Fluxes of CH₄ (a) and N₂O (b) in 2019–2023. The error bars denote standard error. Note the different scale in the y axis in panel (a) for the latter half of the period.

Download

3.3 CH₄ fluxes

Hourly fluxes of CH₄ varied between −50 and 30 µg CH₄-C m⁻² h⁻¹ during the first half of the experimental period (Fig. 3). In conjunction with the rise in the WTD, the values of hourly fluxes increased and varied between −40 and 900 µg m⁻² h⁻¹ during the latter half of the period. The annual flux of CH₄ varied from −1.6 to 11 kg CH₄-C ha⁻¹ yr⁻¹, with an increasing trend towards the end of the measurement period (Table 3; Fig. S3). When the mean annual WTD was below 40 cm, the soil was mainly consuming CH₄, but the consumption tended to change to emissions as the WTD rose (Fig. 2). Variation in the annual cumulative fluxes of each plot was explained by the WTD (p=0.015) but not by crop type (Table S4). The increasing trend between the years 2019 and 2022 was also shown in the mixed-model analysis as the year had a significant effect (p=0.0003) and the effect of WTD decreased with time (years).

3.4 N₂O fluxes

Hourly fluxes of N₂O varied between −3 and 2500 µg N₂O-N m⁻² h⁻¹ during the 4 years, with the highest emissions during the first 4 months (Fig. 3). Annual fluxes of N₂O varied from 1.7 to 33 kg N₂O-N ha⁻¹ yr⁻¹ (Table 3). The emissions declined in time (Fig. S3), especially in the case of set-aside and willow, whereas those of forage did not show such a trend (Table 3). WTD explained the variation in N₂O fluxes well (p=0.015) (Fig. 2; Table S4). There were some interactions between the year and the crop, but the crop type did not affect N₂O emissions systematically between years. Annual N₂O fluxes of the forage and willow treatments differed in the whole time series (p=0.026).

3.5 Global warming potential

The total emissions expressed as CO₂ equivalent ranged from 23 to 39 Mg ha⁻¹ yr⁻¹ with NEE as the CO₂ component and from 27 to 52 Mg with the C export in harvest taken into account in the forage and set-aside treatments (Table 3). For willow, the annual NECB cannot be calculated, but based on the 4-year estimate on carbon binding in the biomass and carbon exported in the harvest divided to a single year, together with the average annual soil respiration and N₂O and CH₄ fluxes, the average annual climate impact of willow cultivation was 10.2 Mg ha⁻¹ yr⁻¹.