Annual balances and extended seasonal modelling of carbon fluxes from a temperate fen cropped to festulolium and tall fescue under two‐cut and three‐cut harvesting regimes

This study reports the annual carbon balance of a drained riparian fen under two‐cut or three‐cut managements of festulolium and tall fescue. CO2 fluxes measured with closed chambers were partitioned into gross primary production (GPP) and ecosystem respiration (ER) for modelling according to environmental factors (light and temperature) and canopy reflectance (ratio vegetation index, RVI). Methodological assessments were made of (i) GPP models with or without temperature functions (Ft) to adjust GPP constraints imposed by low temperature (<10 °C) and (ii) ER models with RVI or GPP parameters as biomass proxies. The sensitivity of the models was also tested on partial datasets including only alternate measurement campaigns and on datasets only from the crop growing period. Use of Ft in GPP models effectively corrected GPP overestimation in cold periods, and this approach was used throughout. Annual fluxes obtained with ER models including RVI or GPP parameters were similar, and also annual GPP and ER fluxes obtained with full and partial datasets were similar. Annual CO2 fluxes and biomass yield were not significantly different in the crop/management combinations although the individual collars (n = 12) showed some variations in GPP (−1818 to −2409 g CO2‐C m−2), ER (1071 to 1738 g CO2‐C m−2), net ecosystem exchange (NEE, −669 to −949 g CO2‐C m−2) and biomass yield (556 to 1044 g CO2‐C m−2). Net ecosystem carbon balance (NECB), as the sum of NEE and biomass carbon export, was only slightly negative to positive in all crop/management combinations. NECBs, interpreted as emission factors, tended to favour the least biomass producing systems as the best management options in relation to climate saving carbon balances. Yet, considering the down‐stream advantages of biomass for fossil fuel replacement, yield‐scaled carbon fluxes are suggested to be given additional considerations for comparison of management options in terms of atmospheric impact.


Introduction
Cultivation of natural peatlands is initiated by drainage, which exposes the peat to aerobic conditions leading to continuous mineralization and CO 2 emissions (Maljanen et al., 2010;Tiemeyer et al., 2016). In addition, the drainage leads to physical peat subsidence resulting from consolidation, shrinkage and compaction (Berglund & Berglund, 2010). In Denmark, the majority of peatlands were drained within the last century for agricultural use, but due to continued mineralization, the total area of agricultural peatlands is decreasing as the soils lose their high carbon (C) content (Greve et al., 2014). To extend the lifespan of peat soils (or their duration in agricultural use), there is a need to pursue new management options to reduce CO 2 emissions from the soil. In this respect, Danish ministerial incentives for extensive management of peatlands were launched in 2015 which call for practices of ceasing tillage, fertilization and pesticide application as well as optionally discontinuing drainage (Ministry of Environment and Food of Denmark, 2015). Such practices may transform the productive cultivated peatlands into extensive grasslands with lower CO 2 emissions. However, another management practice combining productivity and lower CO 2 emissions comprise the cultivating of flooding-tolerant perennial energy crops which do not require annual ploughing and sowing, but still require at least some fertilization and occasional soil tillage (Kandel et al., 2013a). Yet, the climate impact of these inputs may be offset by down-stream greenhouse gas (GHG) savings in terms of replacement of fossil fuels by bioenergy derived from the perennial crops (Shurpali et al., 2010). Indeed, several perennial grasses can be cultivated on rewetted or poorly drained peatlands, which offer a potential to reduce soil C emissions due to attenuated soil respiration (Karki et al., 2016).
Festulolium (9Festulolium) and tall fescue (Festuca arundinacea Schreb.) are suitable perennial grasses to cultivate on peatlands as feedstocks for bioenergy such as biogas production from anaerobic digestion (Sepp€ al€ a et al., 2009;Ambye-Jensen et al., 2013). Festulolium is a hybrid or a hybrid derivative between any species of fescue (Festuca) and ryegrass (Lolium) designed for their combined complementary characteristics (Ghesqui ere et al., 2010;Østrem et al., 2013). These grasses are known for quick establishment, aggressive growth and weed suppression, regrowth capacity, winter survival and high yield. Spring growth of festulolium in particular is earlier than, e.g., perennial ryegrass, and therefore, festulolium can be harvested earlier (Helgad ottir et al., 2014). Both festulolium and tall fescue are native to northern Europe and have large genetic variation which is a useful basis for breeding of cultivars with desired characteristics (Ghesqui ere et al., 2010).
When perennial grasses are cultivated for biogas production, harvesting management is important for optimal quantity and quality of the biomass (Kandel et al., 2013b,c;Wahid et al., 2015). Generally, higher quantity of biomass with better digestibility can be achieved with increased frequency of harvesting, but such managements can increase the environmental impact as frequent harvesting also increases, e.g., the fertilizer demand and fuel consumption related to harvest machinery. To be economically sustainable, the increased cost should at least be offset by the value of additional biomass from frequent harvesting management, but also the environmental cost should be offset by better bioenergy production. Results from studies with reed canary grass (RCG), which is another potential bioenergy crop, have shown that harvesting under two-cut (2C), rather than one-cut management significantly increased biomass yield, CO 2 uptake and total ecosystem respiration (ER) from a drained peatland in Denmark (Kandel et al., 2013a). Likewise, J€ arveoja et al. (2016) reported that fertilizer application in RCG cultivated on peat soil increased biomass production, gross primary production (GPP) and the ratio of autotrophic to heterotrophic respiration. In contrast, only little is known about how improved agronomical managements such as harvest frequency can improve CO 2 uptake and yield of festulolium and tall fescue and thereby improve the GHG balance of these bioenergy crops.
In the present study, we measured the annual balance of CO 2 fluxes in 2C and three-cut (3C) managements of festulolium and tall fescue. The objective was to assess and quantify the effect of harvest frequency on annual CO 2 fluxes, biomass yield and carbon balances. Furthermore, the study also focused on enhancement of empirical modelling of CO 2 fluxes using robust seasonal approaches. This included the use of temperatureadjusted GPP models, to correct for photosynthetic constraints imposed by low temperature, and the use of GPP-derived vegetation proxies in ER models.

Study site and soil properties
The study site was a temperate fen peatland in the Nørre a river valley, Denmark (56°27 0 28 00 N, 9°40 0 33 00 E), which was drained in the early twentieth century for agricultural use. From 1985 to 2009 the site was cropped to rotational grasses and spring barley; from 2009 to 2013 the cropping system included perennial RCG (Kandel et al., 2013d), and in 2013 the field site was ploughed and sown with festulolium (cv. Hykor) and tall fescue (cv. Kora).
Details of the study site are presented elsewhere (Kandel et al., 2013d). In brief, peat depth was more than 1 m with highly decomposed topsoil and less decomposed subsoil. The topsoil (0-20 cm) had a bulk density of 0.29 g cm À3 , total organic C content of 36.3%, total nitrogen (N) content of 3.0% and pH KCl of 6.6. Long-term annual averages  of air temperature, precipitation and global radiation were 7.9°C, 650 mm and 3540 MJ m À2 , respectively.

Experimental design and harvest managements
The study design comprised nine blocks (18 9 24 m) in a 3 9 3 Latin square ( Fig. 1) of which six blocks with festulolium and tall fescue were included in the present study (the three other blocks were cultivated with willow). Festulolium and tall fescue were sown in late August 2013 and fertilized in spring 2014 with standard mineral fertilizers (80:11:38 kg N:P:K ha À1 ). The grasses were harvested in mid-June 2014 followed by mineral fertilization (same rate as in spring), and a second (final) harvest was carried out in late August 2014 (Fig. 2).
The current study was carried out in the second year biomass stands from 01 March 2015 to 29 February 2016 as a splitblock experiment where each block was divided into four plots (18 9 6 m) to impose different harvest managements. The harvest managements included a 3C system with first harvest on 18 May, second harvest on 17 July and third (final) harvest on 30 September. The other three plots in each block were managed as 2C systems differing in the date of the first cut . However, only one 2C system was included in this study, i.e. with first harvest on 15 June and second (final) harvest on 30 September. The 3C and 2C managements received similar rates of fertilization (80:16:60 kg N:P:K ha À1 ) in spring 2015. After each harvest (except final), the harvested areas were supplemented with standard mineral fertilizer at a rate of 80:0:100 kg N:P:K ha À1 (Fig. 2).
At all harvest events, biomass was collected from a central area of 30-45 m 2 in each plot, regarding the area outside the central area as a border. Harvesting was carried out using a sickle-bar harvester (Model GF 3; Grillo SPA, Cesena, Italy). Manual harvesting by hand cutting was carried out for grasses growing in the collars (0.30 m 2 ) installed to measure CO 2 fluxes (described in the following section). All fresh biomass in the collars was weighed and oven-dried at 60°C for determination of dry weight. Non-growing season biomass inside the collars was quantified after the last CO 2 flux measurement on 29 February 2016. Yields were converted to biomass C assuming 45% carbon in dry biomass.

Measurements of CO 2 fluxes
CO 2 fluxes were quantified using a closed chamber approach from 01 March 2015 to 29 February 2016, i.e. covering both the crop growing period (March to September 2015) and the nongrowing period (October 2015 to February 2016). Two weeks prior to the first gas flux measurements, a PVC collar (55 9 55 cm) was inserted to 10 cm depth in each of the 12 plots to measure CO 2 fluxes (Fig. 1). The collars, which remained in the same position throughout the study period, had a 4 cm wide outer flange that was kept parallel to the soil surface to support the top chamber used for flux measurement (Elsgaard et al., 2012). Fluxes were generally measured in 2 weeks intervals, but more frequently during peak growth periods; in total, 29 flux measurement campaigns were conducted during the study year. Fluxes were measured using a closed transparent plexiglass chamber (60 9 60 9 41 cm) equipped with inside and outside temperature sensors, a PAR (photosynthetically active radiation) sensor, a temperature control unit, and cooling and mixing fans (see chamber details in Elsgaard et al., 2012). An extension with same dimension as the top chamber was used when the plant height exceeded the chamber height. The chamber was connected to a portable analysing unit (by 4 mm inlet and outlet tubing) which included a pump, a filter, a portable Li-840 infrared gas analyser (Li-Cor Inc., Lincoln, NE, USA) and a datalogger (CR-850; Campbell Scientific, Logan, UT, USA). Chamber headspace air was circulated through the gas analyser to measure CO 2 and H 2 O concentrations at 1 s intervals.
At each measurement date, four CO 2 fluxes were measured consecutively at each collar at different levels of PAR imposed by shading (Elsgaard et al., 2012). First, CO 2 flux was measured for 1 min under fully transparent conditions; this CO 2 flux was considered as the net ecosystem exchange (NEE) of CO 2 at the prevailing PAR. Immediately after this measurement, the chamber was covered with a netted fabric to block~50% PAR. Plants inside the chamber were allowed to adapt to the new PAR for ca. 1 min before the subsequent 1 min flux measurement. During the adaptation period, one side of the chamber was lifted 15-20 cm from the collar to allow ambient wind to vent the air inside the chamber back to atmospheric CO 2 level. Fig. 1 Layout of the field experiment. Nine blocks (18 9 24 m) in a Latin square were cultivated with festulolium, tall fescue and willow; the present study targeted festulolium and tall fescue. Each block of festulolium and tall fescue was divided into four plots (6 9 18 m) in a split-block design to impose different harvest managements, here including two cuts (2C) or three cuts (3C) in a year. Black squares at the northern end of each 2C and 3C plot indicate the collars (0.55 9 0.55 m) installed for CO 2 gas flux measurements. Fig. 2 Timeline of harvest and fertilizer managements in festulolium and tall fescue plots. Arrows indicate sowing, vertical grey marks indicate fertilizer dates, and vertical black marks indicate harvest dates. Three-cut (3C) and two-cut (2C) harvest managements were imposed in the study year (March 2015 to February 2016; unhatched white area); prior to this date the management was similar in all plots (hatched grey area). See text for details on rates of mineral fertilization.
After the second flux measurement, the procedure was repeated using two layers of netted fabric to block~75% PAR. Finally, the chamber was covered with an opaque cover to block 100% PAR, and the flux was measured after 1 min dark adaptation; this flux measured under dark conditions represented the total ER. The chambers were only partially lifted from collars (15-20 cm from bottom) between the four consecutive flux measurements in order to obtain a steady light response of photosynthesis at decreasing irradiance steps (i.e. no sudden exposure to full PAR) and thereby to avoid a photosynthesis overshoot in the shaded NEE and ER measurements (cf. G€ orres et al., 2014). All fluxes were measured at midday between 10:00 and 13:00. Flux measurement dates were normally fixed in advance and carried out under the prevailing cloud conditions (sunny, partially cloudy or overcast days). However, each individual flux was measured under conditions of nearly constant PAR for 2 min (1 min adaptation and 1 min flux measurement). In partially cloudy days, this was achieved by observing cloud and wind speed just before a flux measurement.
Fluxes were calculated by linear regression using the MAT-

Ratio vegetation index
Growth of grasses was monitored non-destructively by measurement of canopy light reflectance using a SpectroSense 2+ fitted with four SKR1850 sensors (Skye Instruments, Powys, UK) as described by G€ orres et al. (2014). The sensors measured incoming and reflected radiations at red (656 nm) and nearinfrared (778 nm) regions, and were held at nadir angle at a height of 1 m above the soil surface to cover the ca. 0.3 m 2 area inside the collars. Ratio vegetation index (RVI) was calculated as the ratio of reflectance coefficients in the near-infrared and red regions (G€ orres et al., 2014). The RVI was measured prior to all flux measurements and more frequently (2-3 times in a week) between the first and the final harvests to quantify rapid growth of aboveground biomass after the harvest events.

Environmental variables
Volumetric water content (VWC) at 0-20 cm soil depth was measured manually inside each collar on flux measurement days using permanently installed time domain reflectometry (TDR) probes (Thomsen et al., 2007). Soil temperature at 1, 5 and 10 cm depth were recorded continuously (1 h intervals) in the study area using HOBO soil temperature sensors (Onset Computer Corporation, Bourne, MA, USA). Groundwater table (GWT) was measured in a permanently installed piezometer tube; readings were carried out manually from March to mid-May 2015 (at 1-2 weeks intervals) and then continuously during the remaining study period (PTX 1730 pressure transducer sensor; GE Druck, Leicester, UK). Continuous measurements of air temperature and PAR (for extrapolation and upscaling purposes) were obtained from a climate station at Aarhus University in Foulum ca. 7 km from the field site.

Modelling and cumulative CO 2 fluxes
The CO 2 fluxes were partitioned into GPP and ER. Gross primary production was estimated as GPP = NEEÀER, where the sign of ER is positive and sign of GPP is negative; hence, the sign of NEE is positive when there is net flux of CO 2 to the atmosphere. Cumulative annual fluxes were derived from separate modelling of GPP and ER.
For GPP modelling, NEE fluxes measured at each series of full and shaded PAR were corrected for (measured) ER to calculate GPP at the respective PAR levels. GPP was then modelled using a rectangular hyperbolic saturation curve (Thornley & Johnson, 1990) extended to include RVI as a proxy of photosynthetically active biomass (Kandel et al., 2013a): where a (lg CO 2 lmol À1 photon) and GPP max (lg CO 2 m À2 s À1 ), respectively, are the initial slope and asymptote of the light response in the original light response model (Thornley & Johnson, 1990). Parameters were derived by fitting data from all campaigns in one model run as described by Kandel et al. (2013d). The estimated parameters were used to extrapolate annual GPP for the whole study period (01 March 2015 -29 February 2016) using hourly PAR and linearly interpolated RVI measurements at hourly time steps. The GPP model was improved by including a temperature dependence function (Ft) of photosynthesis to account for photosynthetic constraints imposed by low air temperature (e.g. Wu et al., 2015): To scale Ft, the minimum air temperature was set to À2°C as we had observed net photosynthesis down to this temperature during winter and Ft was then assumed to: (i) increase linearly (Ft-Lin) with air temperature to a value of 1 at ≥10°C (Andersen et al., 2013) or (ii) follow an exponential increase (Ft-Exp) between À2°C and 10°C (Fig. 3).
For ER modelling, a temperature modifier function used in the RothC and ECOSSE soil respiration models (Jenkinson et al., 1987;Smith et al., 2010) was applied, but also extended with RVI to model the influence of biomass on ER (Kandel et al., 2013d). Thus, a two-parameter ER model representing basal soil respiration and biomass respiration was developed as: where R b represents a parameter related to basal soil respiration, T air is the air temperature (°C) during chamber enclosure, and b is a scaling parameter for RVI. Other numerals follow the RothC respiration model (Jenkinson et al., 1987). We also fitted the model using soil temperatures measured at different depth (1, 5 and 10 cm depth) but found the best fit using T air , which is reported here.
Total annual ER was calculated following a similar extrapolation approach as for GPP based on hourly air temperatures and linearly interpolated RVI. Annual NEE of CO 2 for each collar was then calculated from the cumulated annual GPP and ER as NEE = GPP + ER.
Modelling uncertainties of annual GPP and ER fluxes were calculated as separate extrapolations of estimated model parameters AE standard errors for each collar (Elsgaard et al., 2012;Kandel et al., 2013d). Modelling uncertainties in annual NEE was then calculated using the quadrature rule of error propagation.

Flux modelling on partial datasets
To test the sensitivity of cumulative CO 2 fluxes to the frequency and timing of measurement campaigns, modelling was also carried out on datasets from alternate measurement campaigns (Alt model), i.e. including only data from every second measurement campaign. Similarly, to test whether a model developed just from growing period data could effectively predict annual fluxes, model parameters were derived from datasets comprising only the growing period (Grow model), i.e. from the measurement campaigns from March to final crop harvest in September 2015. The parameters and annual balances obtained with Alt and Grow models were compared with those obtained by fitting the models to full datasets.

GPP parameters as biomass proxy for ER modelling
The model-derived GPP at PAR of 1000 lmol photons m À2 s À1 (i.e. GPP 1000 ) was tested as a proxy for biomass to be used in ER modelling. Based on the flux measurements at 0, 25, 50 and 100% of ambient PAR, the response of GPP fluxes to light for individual measurement days was modelled for each collar using the rectangular hyperbolic saturation curve of Eqn (1) without RVI. Thereafter, GPP 1000 of biomass in each measurement campaign was calculated. Overall, the light response curves showed a significant relation between GPP and PAR; yet in few occasions (e.g. nongrowing period on overcast days), convergence of the model algorithm failed and linear relations between GPP and PAR were obtained. Under such conditions, only the initial slope parameter [a in Eqn (1) without RVI] was estimated as the saturation parameter (GPP max ) was set to a value close to the previous campaign based on the best model fit (R 2 ).
For use in ER flux extrapolation, GPP 1000 was set to zero on harvest dates and then linearly interpolated at hourly time steps between the GPP 1000 values derived at measurement campaigns. The annual ER obtained using either RVI or GPP 1000 as biomass proxies was then compared. Prior to modelling, both biomass proxies were normalized in relation to the maximum value for individual collars.

Model evaluation
The modelled versus measured GPP and ER fluxes were evaluated using 1:1 scatter plots and statistical indicators such as mean bias, correlation coefficient (r) and paired t-tests (Haefner, 2005). In addition, the Nash-Sutcliffe modelling efficiency (ME) was calculated (Nash & Sutcliffe, 1970); this ME can range from À∞ to 1, where values closer to 1 indicate higher accuracy of the model.

Statistical analysis
Biomass yield at different harvest dates from the same plot was summed to obtain the total annual biomass yield for statistical analysis. Effects of crop type (festulolium, tall fescue), harvesting management (2C, 3C) and their interaction on biomass yield, estimated flux parameters and annual fluxes were tested using analysis of variance (ANOVA) considering block as a random variable. Mean separation was carried out using Fisher's LSD method at 5% level of significance. The statistical analyses were carried out by a mixed model in PROC MIXED (SAS 9.4; SAS Institute, Cary, NC, USA). Biomass yield in plots and collars was compared for each crop/management combination by t-test at 5% level of significance. Results are presented as mean AE SE (standard error) with n = 3, where n is the number of replicated blocks for each crop/management combination. For the individual measured fluxes and supporting variables, ANOVA was not performed as the temporal effects were obvious, and thus, only seasonal trends are presented.
Flux models [Eqns (1)- (3)] were fitted to data from individual collars using the dynamic curve fit function in SigmaPlot 11 (Systat Software Inc., San Jose, CA, USA), which is based on the Levenberg-Marquardt algorithm. The parameters obtained for individual collars were then used for model extrapolations for the particular collar. Pearson's correlation coefficients (r) were presented to indicate the correlation of environmental factors and RVI with measured fluxes.

Climate and groundwater conditions
The mean annual air temperature during the study year (8.4°C) was 0.5°C higher than the long-term mean Fig. 3 Temperature dependence function (Ft) for scaling of gross primary production fluxes. Ft was set to 0 below À2°C and to 1 at and above 10°C; at À2 to 10°C either a linear (Ft-Lin) or an exponential increase (Ft-Exp) response function was assumed. (Fig. 4a). The largest anomaly was December 2015 which was 4.6°C warmer in the study year than normal years. Total global radiation during the study period (3530 MJ m À2 ) was similar to the long-term mean, but April and August received more and May received less radiation than in normal years (Fig. 4b).
The study year was rather wet with total annual precipitation (890 mm) being 240 mm higher than the long-term mean. Notably, May-June and November-December received more precipitation than in normal years (Fig. 4c).
Precipitation showed a fairly uninterrupted pattern during the growing period (Fig. 4d), and accordingly, the grasses did not show any symptoms of water deficit stress. The GWT generally fluctuated at 10-40 cm below the soil surface during the growing period, but was closer to the surface (0-16 cm) in the non-growing period (Fig. 4e) which is normal at the study site.

Environmental and soil conditions during flux measurements
The flux measurements were conducted throughout the study year and thus covered air temperature ranging from 0 to 24°C (Fig. 5a, b). As flux measurements were carried out under prevailing weather conditions (sunny, partially cloudy or overcast days), the ambient PAR also covered a range of intensities, e.g., from 350 to 1820 lmol photons m À2 s À1 during the growing period (Fig. 5c, d). Normally, the mean PAR in a measurement day was similar among the four crop/ management combinations, but differences could be observed when measurements were carried out on partially cloudy days. Anyway, the starting point of flux measurements was randomized in successive campaigns to avoid a bias from diel PAR variation among the treatments. As precipitation was rather evenly distributed throughout the year, VWC was rather stable and close to saturation throughout the study year, but decreased slightly in few occasions in spring and late summer (Fig. 5e, f). VWC was similar across all crop/ management combinations, and there were no obvious effects of harvesting.

Biomass development
Ratio vegetation index of festulolium was slightly higher (0.5-0.8 units) than for tall fescue in the beginning of the growing period from March to mid-April 2015 (Fig. 5g, h). This is consistent with a higher cold tolerance and better winter growth of festulolium. After mid-April, RVI increased more rapidly in tall fescue plots and remained slightly higher than for festulolium. Peaks of RVI were reached (on average) 22 days after harvest events and then declined with similar rates for both crops. After the final harvest in autumn, RVI increased again and fluctuated with temperature during the non-growing period (as seen notably during the warm December). During the non-growing period after final harvest, RVI of festulolium was slightly higher than for tall fescue.

Dynamics of CO 2 exchange
The CO 2 fluxes followed expected seasonal patterns with large fluxes during the crop growing period which decreased during the non-growing period. The midday NEE fluxes of CO 2 were negative except in few occasions just after crop harvest (Fig. 5i, j). Highest net sink strengths (~1600 lg CO 2 m À2 s À1 ) were observed when flux measurements were carried out during sunny days and peak of RVI. The maximum sink strength was similar for the four crop/management combinations. The midday measurements of NEE fluxes were slightly negative during most of the non-growing period as the plots contained regrown biomass after final harvest in autumn and the soil respiration was restricted due to high GWT and low temperature.
Ecosystem respiration also followed a strong seasonal trend (Fig. 5k, l) and correlated well with air temperature (r = 0.87-0.93) and RVI (r = 0.76-0.88), indicating that temperature and green biomass were the major driving factors of ER. The large contribution of aboveground biomass to ER was further evidenced as sharp drops in ER after harvest events. VWC correlated strongly with temperature and RVI and including VWC in the ER model did not significantly improve the model performance (data not shown).
The GPP fluxes showed a clear dependence on weather conditions and plant phenology (Fig. 5m, n). When only GPP measurements taken on sunny days were considered, the fluxes mirrored RVI indicating a strong correlation between crop development and GPP fluxes. However, GPP fluxes were constrained by low PAR on cloudy days and thus large fluctuations were seen in growing period GPP fluxes even for consecutive measurement campaigns within a week. Rates and dynamics of GPP fluxes from both grasses under same harvest management were mostly similar; the few observed differences were mainly related to differences in PAR during flux measurements on partially cloudy measurement days (Fig. 5c, d).

Biomass yield
Total annual biomass yield from the plot area in all crop/management combinations was similar (P > 0.05) with an overall mean of 18.3 Mg dry biomass ha À1 yr À1 (Fig. 6a). Overall, the average biomass yield inside the collars (19.3 Mg dry biomass ha À1 yr À1 ) was slightly higher than that in plots. However, the difference was not consistent in all crop/management combinations, and statistically the yield in plots and collars was similar in all individual crop/management combinations (t-test, P > 0.05). Although regrowth biomass after the final harvest remained green throughout the winter, aboveground biomass yield during the non-growing period (October 2015-February 2016) was modest (Fig. 6b). Festulolium yield during the non-growing period was slightly higher than tall fescue.

CO 2 flux modelling
Effect of temperature functions (Ft) in GPP flux modelling. The GPP models without Ft (No Ft) had excellent performance (ME > 0.92), but still adjustments by temperature functions (Ft-Exp and Ft-Lin) improved model performance for all models run for individual collars (Table S1). Ft-Exp performed slightly better than Ft-Lin, and paired t-tests always showed no differences in measured and modelled GPP fluxes for models run for individual collars when Ft-Exp was used in model fitting (Table S1). Combined 1:1 plots of measured vs. modelled GPP for all crop/ management combinations clearly demonstrated lower bias obtained for fluxes at <10°C when applying a temperature function for GPP adjustment (Fig. 7a-c).
Without temperature function adjustment, the modelling resulted in 12% higher annual GPP compared to the Ft-Exp approach as shown in Fig. 7d (P < 0.05; paired t-tests on estimates for individual collars). Specifically, the GPP model without temperature function estimated about 95% higher GPP in the non-growing period than the GPP model with Ft-Exp (data not shown). This high GPP estimate caused a low ratio of biomass-to-GPP (biomass/GPP) in the non-growing period which averaged only to 0.15; using the Ft-Expadjusted GPP, the biomass/GPP ratio increased to 0.28 (Fig. 7e). As the overestimation of non-growing season photosynthesis was corrected by temperature adjustments, the GPP model with Ft-Exp was used for further comparisons of crops and managements.
Resulting GPP and ER model performance. The resulting GPP model (i.e. with Ft-Exp temperature adjustment) run for individual collars had excellent calibration fit (Fig. 8a, b) as signified by high MEs of 0.94-0.97 (Table S1). The model parameters (GPP max and a) were always significant (P < 0.05), and both parameters were statistically similar among the four crop/harvest combinations, i.e., showing less than 12% deviation between individual collars (Table S1).
The ER model also showed a very good fit to the measured data (Fig. 8c, d) with significant model parameters (P < 0.05) and MEs ranging from 0.85 to 0.95 for individual collars (Table S2). Paired t-tests showed no significant differences in measured and modelled ER fluxes for all models run for individual collars. As the flux models included biomass as a major controlling factor, the effects of rapid changes in biomass status (e.g. by harvest events) can be seen as lowered GPP fluxes and thereby decreased net uptake (NEE) of CO 2 (Fig. S1).

Comparison of annual balances
Festulolium under 3C management had the highest GPP and ER and tall fescue under 2C management had the lowest GPP and ER (Fig. 9a), but differences (and interactions) were not significant (P > 0.05). NEE was also statistically similar among all crop/management combinations and ranged between À750 and À840 g CO 2 -C m À2 yr À1 . Integrating the NEE of CO 2 -C and the export of harvested biomass carbon to represent the net ecosystem carbon balance (NECB) showed that the management systems acted as weak sources of carbon except for the tall fescue 2C system that was a weak sink of carbon. However, as for biomass yield and NEE, NECB was also statistically similar among the four crop/management combinations (Fig. 9b).
The spatial variation (uncertainty) at block scale in C fluxes by GPP, ER and biomass export was small, but nevertheless translated into relatively high uncertainties in the estimates of notably NECB (Fig. 9b). Likewise, the magnitude of uncertainty of GPP and ER fluxes related to the modelling approach was low (as the model performance was high), but the propagation of uncertainties in GPP and ER resulted in a relatively higher uncertainty of NEE fluxes (Fig. 9a).
Flux modelling with partial datasets. The model parameters obtained for modelling on full datasets and Alt datasets (i.e. including only every second campaign) were similar for both GPP and ER models (Tables S1  and S2). Likewise, although the number of measurements in Alt datasets was only half of the full datasets, the Alt datasets covered a similar range of fluxes and driving variables as the full datasets. The mean MEs of fluxes derived from the Alt datasets were 0.97 for GPP and 0.91 for ER; with the full dataset, the corresponding MEs were 0.97 and 0.90, respectively (Tables S1 and S2).
The GPP parameters estimated for the Grow (i.e. growing period) datasets were also similar to parameters for the full datasets. The Grow datasets covered flux measurements under large PAR ranges obtained by shading and different cloud conditions. Likewise, the Grow datasets covered large RVI ranges provided by harvest and regrowth during the growing period. MEs of GPP models developed with Grow datasets (average 0.94) were as strong as MEs for models developed with the full datasets (average 0.95) (Table S1). However, the performance of ER models based on the Grow datasets (average ME = 0.77) was weaker than for models with the full datasets (average ME = 0.90) (Table S2). This was attributed to a narrower temperature range in the Grow datasets (13-24°C) than in the full datasets (0-24°C).
The annual fluxes of GPP obtained with Alt and Grow datasets deviated by <5% and 1%, respectively, from the results for the full datasets, and these differences were not statistically significant (Fig. 10a, b). Annual ER fluxes obtained with full and partial datasets also were similar (Fig. 10c, d). Similarly, differences between measured and modelled fluxes were not significant for all Alt models developed for individual collars (P > 0.05; paired t-tests, Table S2).
GPP 1000 as biomass proxy for ER modelling. The GPP 1000 and RVI (both normalized relative to their maximum values) showed strong intercorrelations for measurements from all collars, indicating that GPP 1000 may potentially substitute RVI as a proxy of standing green biomass (Fig. 11a). Indeed, substituting GPP 1000 for RVI in the ER models resulted in very good model fits (average ME = 0.90) similar to the models using RVI (Table S1). Also the ER model parameters obtained with GPP 1000 and RVI were similar (Table S2). Thus, the annual ER estimated with GPP 1000 as biomass proxy was similar to annual ER estimated with RVI as biomass proxy (Fig. 11b).

Biomass development and yield
To our knowledge, the biomass yield observed in this study (18.3-19.3 Mg ha À1 ) represents the highest yield of festulolium and tall fescue so far reported in northern (e, f) volumetric water content (VWC) measured at 0-20 cm soil depth and normalized in relation to the maximum during flooded winter conditions; (g, h) ratio vegetation index (RVI) as a proxy of green biomass. RVI was measured more frequently (2-3 times in a week) between the harvest events, but for clarity only measurements taken at the same time as CO 2 flux measurements are shown. Lower panels show the measured CO 2 fluxes of (i, j) net ecosystem exchange (NEE) at prevailing PAR (measurements without shrouds), (k, l) ecosystem respiration (ER) and (m, n) gross primary production (GPP) calculated as NEEÀER. CO 2 fluxes were measured at midday from 10:00 to 13:00 at prevailing PAR. The dashed vertical lines indicate dates of grass harvest. Error bars represent the spatial variations at block scale (standard error, n = 3). Unidirectional error bars are shown for clarity.
Europe (Sepp€ al€ a et al., 2009;Ambye-Jensen et al., 2013Østrem et al., 2013Butkut_ e et al., 2014). However, similar yields were observed in the same growing period in a field trial on mineral soil close to our study site (Manevski et al., 2016). In the current study, both festulolium and tall fescue stayed green even under flooded and frozen soil conditions during winter indicating their potential as biomass crops under harsh environmental conditions.
In flux studies with closed chambers, plants inside the permanent collars can be disturbed by repeated chamber deployment and consequently biomass yield can be negatively affected (Maljanen et al., 2001). However, careful handling and hardiness of the cultivated grasses prevented such effects in the current study as slightly, though not significantly higher yields (ca. 1 Mg ha À1 ) were recorded inside the collars than in the surrounding plots. Likewise, dynamics of RVI inside the collars and the plots were similar (data not shown). Thus, the vegetation status in the collars was considered to be representative of the plots in terms of biomass growth and yield, which was a major controlling factor of CO 2 fluxes.

Flux modelling
The GPP model including RVI as biomass proxy was previously presented by Kandel et al. (2013a) and successfully adapted in other studies (J€ arveoja et al., 2016;Karki et al., 2016). In the current study, the GPP model was extended with a temperature function (Ft) in analogy to the commonly used approach to describe the effect of low temperatures on radiation use efficiency of C 3 crops (Andersen et al., 2013). Although the overall fit was satisfactory also for GPP models without Ft, these models were biased especially in winter when photosynthesis was likely to be constrained by low temperature (Fig. 7a). This high GPP estimate caused an unlikely low ratio of biomass/GPP in the non-growing period which averaged only to 0.15 (Fig. 7e). With the Ft-Exp-adjusted GPP, the average biomass/GPP ratio was 0.28 indicating a much more realistic GPP estimate (Trischler et al., 2014) and thus also indicating an overstimation of wintertime GPP by the model without Ft. Including Ft in the current GPP model provided a more robust model, and GPP overestimation in colder periods was effectively corrected.
Greenness indices, such as RVI, can be useful vegetation proxies for ER modelling (Gilmanov et al., 2005;Kandel et al., 2013a,d), but canopy reflectance or other crop development parameters are not routinely included in flux study protocols (e.g. Hoffmann et al., 2015). However, ER models based only on temperature often show a lack of fit to seasonal data due to confounding contributions from biomass dynamics and therefore some seasonal studies include cluster-wise ER modelling for distinct plant phenological stages to overcome such confounding effects (Beetz et al., 2013;Poyda et al., 2016). Another approach has been to use instantaneous GPP measured under ambient PAR as a biomass proxy (Larsen et al., 2007;Kopittke et al., 2013). However, GPP is strongly sensitive to PAR, and coupling ER with GPP therefore also makes ER sensitive to diel PAR fluctuations which may not be generally valid. Yet, a stable photosynthesis measure, which is not affected by diel variations in PAR, can be a useful biomass proxy in ER models, particularly in those studies without other quantification of crop development. As shown in the present study, GPP 1000 , which is a measure of ecosystem photosynthetic capacity, can be such a useful alternative proxy of vegetation status for ER modelling as it does not couple ER to diel PAR changes.
An often stated weakness of manual chamber methods for annual flux estimation is the large data gaps resulting from low sampling frequency as compared to the long (often annual) period of interest, thus making extrapolation crucial. In the current study, for example, the total chamber enclosure time per collar was ca. 5 h corresponding to <0.1% of the full year, even though 29 measurement campaigns were achieved. Therefore, robust models and significant parameters for gap-filling are central for the validity of annual fluxes extrapolated from manual chamber measurements. In the current study, biomass and temperature explained most of the variations in ER fluxes. With the current seasonal models, robust ER and GPP parameters were obtained even though every other measurement campaigns were excluded in the model calibration (Alt datasets); this indicated that it may be more important to capture the full amplitude of driving variables rather than achieving a high frequency of measurements. This was supported by the deviating ER results from modelling with the Grow datasets which represented a higher frequency of measurements, but a narrower range of driving variables (temperature) than the Alt datasets. On the other hand, GPP modelling with the Grow datasets gave similar results as for the Alt datasets as both datasets captured a similar range of the GPP driving variables (PAR and RVI).
Previous field studies have indicated a large range of apparent temperature sensitivities of ER on seasonal scales, which may be related to combined effects of temperature and other driving variables co-varying with temperature (Curiel Yuste et al., 2004;Reichstein et al., 2005). Similarly, ER modelling using temperature from different soil depths may also cause apparently different temperature sensitivities due to time lags between CO 2 production at different soil depths and the CO 2 emission from the soil surface (Phillips et al., 2011;Zhang et al., 2015). In the current study, extrapolation to annual ER based on temperature sensitivity obtained with seasonal data was avoided by using a generic temperature sensitivity function (Jenkinson et al., 1987) as applied in some mechanistic models including RothC and ECOSSE (Smith et al., 2010) as an empirical temperature modifier function.

Annual CO 2 balances
The similar annual estimates of GPP fluxes in all crop/ management combinations aligned with the similar biomass yields. Compared with previous RCG cultivation at the study site under one-cut regimes (Kandel et al., 2013d), average aboveground biomass yields of festulolium and tall fescue were~40% higher although annual GPP was only~20% higher. This was potentially related to higher root allocation of photosynthates in the one-cut RCG system as perennial grasses allocate more photosynthates to roots in senescence stages when shoot growth stops (Warembourg & Estelrich, 2001;Ge et al., 2012;Zhang et al., 2014). However, compared with RCG cultivation under 2C regimes with fertilization after the first harvesting (Kandel et al., 2013a), the current yields of festulolium and tall fescue were at a similar level.
Previous studies of CO 2 fluxes at the same study site also included ER fluxes measured for a complete year with RCG cultivation under one-cut management (Kandel et al., 2013d;Karki et al., 2015). Although the differences in crop type and managements can have large impacts on biomass yields and thereby ER fluxes, comparisons of RCG and the current grasses provide some indication of interannual variability in ER fluxes influenced by climate differences. Annual ER estimates for festulolium and tall fescue (1400 g CO 2 -C m À2 ) were lower than previously reported for RCG (ca. 1900 g CO 2 -C m À2 ) under one-cut management at the study site (Kandel et al., 2013d;Karki et al., 2015), even though the biomass yields of festulolium and tall fescue were higher. The ER difference partly originated from the modelling approach as the modified van't Hoff type model used in the previous studies also estimated an annual ER of 1900 g CO 2 -C m À2 for the current dataset, which is ca. 34% higher than the estimates with the model in this study (data not shown). However, ER in the current study was still lower than in the RCG studies as calculated per unit of biomass yield. Possible reasons for the lower yield-scaled ER in the current study could be related to differences in allocation of photosynthates to above-or belowground parts as the current grasses (under 2C and 3C management) expectedly had lower root-to-shoot ratios than the RCG under one-cut management which remained in a senescence stage for longer time during the growing period (Warembourg & Estelrich, 2001;Zhang et al., 2014); thus, ER from root compartments could be higher in RCG. A lower ER in the current study could also be related to lower temperatures in the period between May to July which may have contributed to lower soil respiration. Thus, although the annual temperature in the current study was higher than in normal years, this was due to high temperature mainly in December when GWT was close to the soil surface and soil respiration thereby was constrained. Indeed, a likely reason for the lower ER in the current study was the difference in GWT which was ca. 15 cm below soil surface (annual average) as compared to an annual average of ca. 40 cm below the surface in the previous RCG study year (Karki et al., 2015). GWT was not raised intentionally in the current study year, but the higher GWT was at least partly related to a poor cleaning of the drainage ditch in the study area in recent years in combination with the high amount of precipitation in the study year. Fig. 8 Comparison of modelled and measured gross primary production (GPP) and ecosystem respiration (ER) for all 12 collars, i.e. representing both (a, c) three-cut (3C) and (b, d) two-cut (2C) managements of festulolium (black circles) and tall fescue (grey circles). The diagonal lines represent 1:1 agreement.

Use of NECB as emission factor to compare crop management systems
Net ecosystem carbon balances were derived by summing the flux-based C balances (NEE of CO 2 -C) and the removal of carbon in harvested biomass (Elsgaard et al., 2012). Although biomass was statistically similar among crop/management combinations, there was a quite large yield range (556-1044 g C m À2 ) among the individual 12 collars (analogous to the plots) and there was a clear tendency of higher biomass yield under 3C managements (Fig. 12). The large differences allowed developing a correlation between accumulated aboveground biomass yield and the cumulative fluxes of carbon from the ecosystems. Based on measurements from individual collars, there was strong positive correlation (r = 0.85) between biomass yield and NECB (Fig. 12a), whereas there was no correlation (r = 0.02) between biomass yield and NEE (Fig. 12b). The correlation between biomass and NECB indicated that the collars with the lowest biomass production were net sinks of carbon whereas collars with the highest biomass production were net sources of carbon (Fig. 12a). The difference was probably related to differences in root-to-shoot ratio as grasses under 2C management remained in senescing stages for a longer period than the grasses under 3C managements, and the proportion of photosynthates allocated to roots in the senescing period possibly increased (Warembourg & Estelrich, 2001;Ge et al., 2012;Zhang et al., 2014). The 3C management, however, allowed the grasses to develop new shoots, and thus, photosynthates allocation to shoots was expectedly enhanced. This assumption was supported by the higher ratio of photosynthates fixed in aboveground parts from the more productive management systems (Fig. 12c).
The higher NECB from more productive management systems align with previous studies, indicating that management systems which suffered from nutrient deficiency or flooding stress had comparably lower NECB compared to managements without those stresses (Kandel et al., 2013a;Karki et al., 2016). Similar positive correlation between NECB and biomass yield was observed in another study which included measurements in eight peatland sites in Denmark cropped with annual and perennial crops (Elsgaard et al., 2012). Thus, when emission factors (EFs) from different cropping systems are compared in terms of yield-corrected flux balances, the possible advantages of higher yields are to some extent concealed. Including fossil fuel displacement by harvested biomass in EF calculations would mitigate this bias towards favouring low-yielding cropping systems; yet such a detailed life cycle assessment was not within the scope of this study. Alternatively, we suggest that management systems can be reasonably evaluated and compared in terms of ER carbon fluxes per unit of yield. Such yield-scaled ER plotted against yields (Fig. 12d) favoured the high-yielding ecosystem in terms of lower relative ER for higher yields. As soil respiration is an important component of ER, using yield-scaled ER for ecosystem comparisons considers the carbon emissions from soils, and can be a balanced approach between economic and environmental concerns whereas comparing management systems in terms of NECB may partly select against improvement of biomass yield.

Role of other emissions in NECB and GHG balance assessment
Previous studies on drained peatlands in Denmark suggested NEE of CO 2 and biomass removal as the major fluxes contributing to the carbon and GHG balance of drained peatlands (Elsgaard et al., 2012;Karki et al., Fig. 9 (a) Annual fluxes of gross primary production (GPP), ecosystem respiration (ER) and net ecosystem exchange (NEE) of CO 2 . (b) Net ecosystem carbon balance (NECB) calculated as the sum of NEE and biomass carbon exported from the ecosystems at harvest events. Thin error bars (unidirectional) represent the spatial variation at block scale (standard error, n = 3), and thick error bars without cap (in panel a) represent the uncertainty of fluxes related to the modelling approach. 2015). Thus, methane (CH 4 ) fluxes were not measured as part of this study since no significant CH 4 fluxes were observed at the study site in a previous year (Karki et al., 2015). However, CH 4 contributions could possibly come from drainage ditches (IPCC, 2014).
According to the default emission factor for drainage ditches in shallow drained peatlands used as grassland in temperate Europe, the study site would emit ca. 2 g CH 4 -C m À2 yr À1 on landscape scale (IPCC, 2014). Similarly, waterborne losses as dissolved organic carbon Fig. 10 Comparisons of (a, b) annual gross primary production (GPP) and (c, d) ecosystem respiration (ER) obtained by modelling and extrapolation of datasets including all measurements (All) versus partial datasets including alternate (i.e. every second) measurement campaigns (Alt), or only measurements (Grow) from the growing period for each combination of festulolium (F) and tall fescue (TF) and two-cut (2C) or three-cut (3C) management. Slope parameters of the linear regression and goodness-of-fits (R 2 ) are shown. The diagonal lines represent 1 : 1 agreement. Fig. 11 (a) Correlation between RVI (ratio vegetation index) and GPP 1000 (GPP at a photosynthetically active radiation of 1000 lmol photon m À2 s À1 ). Both variables were normalized in relation to their maximum. (b) Correlation between annual ecosystem respiration (ER) obtained by modelling and extrapolation using RVI (ER-RVI) and GPP 1000 (ER-GPP 1000 ) as biomass proxies. Data are shown for all 12 individual collars in the four different crop/management combinations, i.e. festulolium and tall fescue under two-cut (2C) and three-cut (3C) management. The diagonal lines represent 1 : 1 agreement.
(DOC) can also contribute to the total C balances (IPCC, 2014). For a tentative assessment, the average DOC emissions from drained peatlands (21 g C m À2 yr À1 ), presented in IPCC (2014), was used in this study. Thus, the estimated contribution of CH 4 and DOC fluxes would only make a modest contribution (23 g C m À2 yr À1 ) to the total carbon balance as compared to NEE of CO 2 (average, À776 g C m À2 yr À1 ) and biomass removal in harvest (average, 868 g C m À2 yr À1 ).
This study primarily focused on measurement and modelling of CO 2 fluxes. However, in terms of GHG balances, the study site is a possible source of nitrous oxide (N 2 O) fluxes due to soils with low C/N ratio and high N content (Klemedtsson et al., 2005). N 2 O fluxes were not measured, but were low during the study year in adjacent plots studied under flooded and natural conditions (T.P. Kandel, unpublished results). These adjacent study plots were dominated by RCG and Poa sp; but crop managements, such as frequency and timing of harvesting and fertilization, were similar to the current 2C plots and biomass yields were also comparable to the current levels. Annual N 2 O fluxes in the adjacent plots under natural conditions were only ca. 0.5 g N 2 O m À2 yr À1 , and significant emissions were recorded only around the time of fertilization in spring and summer (T.P. Kandel, unpublished results). Assuming similar dynamics and magnitude of N 2 O fluxes in the current study (where the 2C fertilizer amount, timing of application and biomass yields were comparable), it was estimated that N 2 O fluxes under 2C and 3C managements would be equivalent to ca. 0.5 g and 0.75 g N 2 O m À2 yr À1 , respectively (i.e. 1.5 times higher in 3C managements due to additional fertilizer). These N 2 O fluxes corresponded to 40 and 60 g CO 2 -C m À2 yr À1 considering 298 times higher global warming potential of N 2 O in comparison with CO 2 (IPCC, 2007). These results therefore indicate that CO 2 fluxes dominated the total GHG balance at this study site, supporting the results of previous studies (Elsgaard et al., 2012;Karki et al., 2015). (c) proportion of cumulated gross primary production (GPP) allocated to aboveground parts (biomass/GPP ratio); and (d) proportion of ecosystem respiration (ER) relative to aboveground biomass yield (ER/biomass ratio). The fitted linear relations represent data for all four crop/management combinations, i.e. festulolium and tall fescue under two-cut (2C) and three-cut (3C) management. Parameters of the linear regressions and correlation coefficients (r) are also shown.
internship students Suman Thapa and Ben McCarthy for their assistance in field works. Finally, we would also like to thank the anonymous journal reviewers for careful reviews of the manuscript and useful suggestions.