Carbon dioxide ﬂuxes and carbon balance of an agricultural grassland in southern Finland

. A signiﬁcant proportion of the global carbon emissions to the atmosphere originates from agriculture. Therefore, continuous long-term monitoring of CO 2 ﬂuxes is essential to understand the carbon dynamics and balances of different agricultural sites. Here we present results from a new eddy covariance ﬂux measurement site located in southern Finland. We measured CO 2 and H 2 O ﬂuxes at this agricultural grassland site for two years from May 2018 to May 2020. Especially the ﬁrst summer experienced prolonged dry periods, which affected the CO 2 ﬂuxes, and substantially larger ﬂuxes were observed 5 in the second summer. During the dry summer, leaf area index (LAI) was notably lower than in the second summer. Water use efﬁciency increased with LAI in a similar manner in both years, but photosynthetic capacity per leaf area was lower during the dry summer. The annual carbon balance was calculated based on the CO 2 ﬂuxes and management measures, which included input of carbon as organic fertilisers and output as yield. The carbon balance of the ﬁeld was –50 ± 68 g C m − 2 yr − 1 and –118 ± 24 g C m − 2 yr − 1 in the ﬁrst and second study year, respectively. We estimated that on average the grassland exceeded the 10 global "4 per 1000" goal to increase the soil carbon content.


Introduction
Conventional and intensive agricultural practices cause significant carbon emissions while diminishing the soil organic matter (SOM) content. This leads to a reduction of soil quality and health (e.g. Houghton and Nassikas, 2017;Le Quéré et al., 2009Lal, 2016;Paustian et al., 2000;Smith, 2008). Currently, agriculture is responsible for more than 10% of the global 15 anthropogenic greenhouse gas (GHG) emissions to the atmosphere (Le Quéré et al., 2017). Soil type and properties, vegetation, climate and weather conditions as well as management practices all have a considerable effect on the carbon fluxes and balances of agroecosystems (Bolinder et al., 2010;Gomez-Casanovas et al., 2012;Jensen et al., 2017;Lorenz and Lal, 2018;Singh et al., 2018). Frequent ploughing, monocropping and intensive use of agrochemicals are the main contributors to the loss of SOM and the resulting carbon dioxide (CO 2 ) emissions from land use Reinsch et al., 2018;Yang et al., 2019). 20 A change from conventional and intensive agricultural practices to regenerative and holistic farm management provides a substantial climate change mitigation potential (Lal, 2016). Increasing the amount of SOM in agroecosystems by applying enhanced management practices, such as lighter tillage, continuous plant cover, rotational grazing, agroforestry, increased The eddy covariance (EC) method is widely used for measuring CO 2 and energy fluxes in different ecosystems and climatic conditions . The high-frequency measurements provided by EC allow a direct quantification and analysis 60 of gas exchange between the ecosystem and atmosphere. The carbon balance calculated from EC data, combined with the additional carbon fluxes caused by management, serves as an important measure for determining the climatic impact of agricultural ecosystems (e.g. Baldocchi, 2003;Baldocchi et al., 2018). However, continuous GHG flux measurements on agricultural sites, especially on mineral soils and grasslands, are still scarce in the Nordic countries (Shurpali et al., 2009;Lind et al., 2020;Jensen et al., 2017). 65 The aim of this study is to determine the magnitude and seasonal dynamics of the carbon balance of a managed forage grassland in southern Finland. In particular, we had three specific research questions: 1. How does the CO 2 exchange and carbon balances vary between the study years? 2. Does the grass photosynthesis indicate occasional drought-related responses? 3. How does the possible carbon sink relate to carbon sequestration objective of "4 per 1000" initiative? 70 For the purposes of this study, we collected field data on the net exchange of CO 2 and H 2 O, soil and vegetation properties and meteorological variables on an agricultural grassland in southern Finland during two years, from May 2018 to May 2020.

Site description
The flux measurements were conducted at the Qvidja farm in southern Finland (60.29550 • N, 22.39281 • E; elevation 5 m) from 75 May 2018 to May 2020 (Fig. 1). The site belongs to the hemiboreal climate zone. From 1981 to 2010, the mean annual air temperature and precipitation at the Kaarina Yltöinen weather station, located 13 km northeast of Qvidja, were 5.4 • C and 679 mm, respectively (Pirinen et al., 2012). The experimental field in Qvidja has mineral soil (clay loam) and it covers 16.25 ha. It was cultivated as forage grassland during the study years. From 2008 to 2016, the field was managed intensively with conventional practices, and it was in annual crop rotation. In 2017, the field management practices were converted towards 80 more sustainable and environmentally friendly farming by increasing the use of organic fertilisers and perennials, restricting the use of pesticides and increasing plant species biodiversity. The current grass was sown as an undergrown species with broad bean in spring 2017. The predominant grass species were timothy (Phleum pratense), meadow fescue (Festuca pratensis) and white clover (Trifolium repens).
Grass was harvested for silage for the first time on 12 June 2018. As the grass cover was fairly sparse later in the summer due 85 to drought, repair seeding was done on 3 September 2018 to restore the drought-induced damage. The seed mixture included 35% of timothy, 30% of rye-grasses (Lolium spp.), 20% of common meadow-grass (Poa pratensis) and 15% of red fescue (Festuca rubra). Timothy, meadow fescue and clover remained as the predominant species also in 2019 and early 2020. On 21 August 2018, the grass was cut at approximately 15 cm, but the yield was left in the field. The second harvest of 2018 occurred Figure 1. Experimental field with the sectors representing the target area. the area covers 3.9 ha. Eddy covariance tower is located in the centre of the sectors. Wind directions from 30 to 140 • were filtered out due to another experimental plot locating in that part of the field. (Orthophoto from National Land Survey of Finland) on 23 September. In 2019, the grass was harvested on 11 June and 20 August. In June 2018, a conventional cutting height of 6 90 cm was used, whereas in the other harvests the grass was cut at 15 cm.
In 2018, the field was fertilised twice, on 16 July and 24 August, with 2800 kg ha −1 and 1800 kg ha −1 of NK-molasses, respectively (Table 1). NK-molasses was a byproduct of the sugar industry. It contained 67% of organic matter (OM) and 4.4% of nitrogen and had the C:N ratio of 9. According to the product information, the molasses included 205 g kg −1 of organic carbon. In addition, it contained potassium and small proportions of sulphur, magnesium, calcium and sodium.

95
In May 2019, the field was fertilised with a mixture of side products from industries of starch potato processing, biowaste processing and ethanol production out of sawdust. This fertilisation mixture contained 65% of OM, 1.3% of nitrogen, 0.2% of phosphorus, 3% of potassium and 0.4% of sulphur, as well as small amounts of calcium, magnesium, zinc, copper and manganese. Approximately 4600 kg ha −1 was applied on the field on 8 May (Table 1). On 26 June after the first harvest, 220 kg ha −1 of mineral fertiliser was applied. This fertiliser contained 23% of nitrogen, 10% of phosphorus and 8% of potassium. which detects the CO 2 and H 2 O mixing ratios, and a three-dimensional sonic anemometer (uSonic-3 Scientific, METEK GmbH, Elmshorn, Germany) to measure wind speed and air temperature. The data were recorded at 10-Hz frequency. The measurement height was 2.3 m. The flow rate was about 12 l min −1 , and the length of the 4-mm stainless steel inlet tube with 2 µm Swagelok sinter was 0.8 m. The gas analyser was calibrated with a zero CO 2 concentration air as a reference gas in May 2018 and March 2020. The micrometeorological sign convention is used throughout the paper, with a negative value indicating 110 the flux from the atmosphere to the ecosystem (net uptake) and a positive value indicating the flux from the ecosystem to the atmosphere (net emission).
Auxiliary meteorological measurements were conducted next to the flux tower. These included soil moisture observations at the depth of 0.1 m (ML3 ThetaProbe sensor, Delta-T Devices Ltd., Cambridge, UK) and soil temperature profile at the depths of 5, 10 and 30 cm (Pt100 IKES sensors, Nokeval Oy, Nokia, Finland). The soil temperature data were collected (SNAP) software. The cloudy, cloud-shadowed and snowy data were filtered out using the scene classification band available 125 in the L2A products.

Eddy covariance data processing
The turbulent fluxes were determined as the covariance between the variations of vertical wind component and gas mixing ratio recorded at 10 Hz. They were calculated as 30-min block averages applying standard procedures, including double coordinate rotation and lag determination based on cross-correlation analysis (Rebmann et al., 2012). The systematic flux loss due to 130 the incomplete frequency response of the measurement system was corrected according to the empirical method described by Laurila et al. (2005).
The EC data from 5 January to 28 March 2019 were affected by technical issues with an inlet filter, which resulted in an erroneous reading of the internal analyser pressure. For this period, the 10-Hz mixing ratios were recalculated from the recorded absorptance data using the instrument-specific calibration functions. The mean CO 2 mixing ratio was set to 410 ppm in these 135 calculations. The following acceptance criteria were applied to screen the 30-min averaged CO 2 flux data: number of spikes in the raw data < 150 of 18,000, relative stationarity of CO 2 flux (Foken et al., 2012) < 50%, mean CO 2 mixing ratio > 380 ppm, variance of CO 2 mixing ratio < 15 ppm 2 between April and September and < 5 ppm 2 between October and March, and wind direction within 0-30 • or 140-360 • . Furthermore, the data were discarded during the periods of weak turbulence and when the flux footprint was not sufficiently representative of the target grassland, as estimated with the footprint model of Kormann and 140 Meixner (2001). For these, we applied a friction velocity limit of 0.06 m s −1 and a cumulative footprint limit of 0.7. The further screening applied to H 2 O fluxes included: H 2 O flux > 0, relative stationarity of H 2 O flux < 50% and variance of H 2 O mixing ratio < 1 (mmol mol −1 ) 2 . After applying these filtering criteria, the coverage of CO 2 and H 2 O flux data accepted for further analysis was 44% and 30% of all the 30-min periods during the two measurement years, respectively. Most of the accepted CO 2 and H 2 O flux data were collected when the wind direction was in the south-southwest sector (Fig. 2).

Soil temperature model
The soil temperature sensor at the depth of 5 cm malfunctioned during the first measurement year, and these data were replaced with values derived from air temperature using the model presented by Rankinen et al. (2004). This model also takes into account the effect of possible snow cover on soil temperature. The following equation was used to obtain 30-min soil temperatures at 5 cm from 8 May 2018 to 3 May 2019: where T t+1 z is the soil temperature at the depth of Z s on the following day, T t z is the soil temperature of the current day, ∆t is the length of the timestep, K T is soil thermal conductivity, which was set to 1 W m −1 K −1 , C A is the apparent heat capacity which is the sum of specific heat capacity of the soil C s = 0.5×10 −6 J m −3 K −1 and specific heat capacity due to freezing and thawing C ice = 4×10 −6 J m −3 K −1 , and T t air is the measured air temperature. The impact of snow cover was taken into 155 account in the last term of the equation where f s is an empirical snow parameter, which was set to 10 m −1 , and D s is the measured snow depth. The model predictions were compared to measurements at the experimental field between June 2019 and May 2020. During summertime, the changes in soil temperature were fairly well captured by the model, whereas in the wintertime, the model tended to create larger changes in temperature than the actual measurements showed (Fig. 3).

160
To calculate CO 2 balances and to conduct further seasonal analysis of the flux components the measured CO 2 flux data (i.e. net ecosystem exchange, NEE) were partitioned to GPP and total ecosystem respiration (R eco ) and gap-filled based on this partitioning: The gap-filled GPP and R eco were calculated with empirical response functions by first fitting these functions to the flux data.

165
R eco was expressed as a function of temperature (Lloyd and Taylor, 1994): where R 0 is the respiration rate (mg m −2 s −1 ) at the reference soil temperature of 283.15 K, T 0 = 227.13 K, T 1 = 56.02 K, and E 0 = 308 K is the long-term ecosystem sensitivity coefficient (Lloyd and Taylor, 1994) that describes the temperature response of soil respiration, and T s is the soil temperature at the depth of 5 cm. GPP was modelled as a function of photosynthetically active radiation (PAR, µmol m −2 s −1 ) and daily effective phytomass index (PI) as: where PI is an empirically-determined variable introduced to describe the seasonal changes in the photosynthetically active 175 vegetation (Aurela et al., 2001), α is the apparent quantum yield (mg µmol −1 ), and GP max denotes the asymptotic CO 2 uptake rate in optimal light conditions (mg m −2 s −1 ). Further details on the PI determination and gap-filling procedure are provided in the Appendix A and B, respectively. Energy fluxes were gap-filled following the description in the Appendix C.
To study the differences in photosynthetic capacity of the grass field between the two growing seasons, daily GP 1200 values were calculated with the estimated α and GP max values, i.e. GPP was normalised to PAR = 1200 µmol m −2 s −1 .

Carbon balance
The carbon balance of the agricultural ecosystem was calculated by adding up the 30-min NEE fluxes, the imported carbon in the form of organic fertilisers and the removal of carbon as harvested biomass: where C H is the amount of carbon in harvested biomass, C F is the amount of carbon in imported fertilisation and m is the 185 total number of timesteps in the period for which the balance was calculated. Biomass was converted to carbon by multiplying the dry weight by 0.42 (Lohila et al., 2004). The following sign convention was used: the carbon imported into the ecosystem corresponds to a negative flux and the carbon removed from the system corresponds to a positive flux.

Uncertainty analysis
The CO 2 balance, which is calculated based on the EC measurements, includes multiple potential error sources. Uncertainties 190 are associated, for example, with the stochastic nature of turbulence and incomplete sampling of large eddies, the performance of instruments and the flux variation caused by the limited area of the target ecosystem . Some of these errors were compensated for in the data processing and screening. Here we included in the uncertainty estimate the most relevant random error sources, i.e. the statistical measurement error (E meas ) and the error caused by gap-filling (E gap ) Aurela et al. (2002): where NEE meas is the filtered 30-min flux, NEE mod is the corresponding modelled NEE (Eqs. 2-4), and n is the number of measured data.
where E GP P and E Reco are the errors of modelled GPP and R eco , respectively. N is the number of gaps in the data.

200
The standard error propagation principle was used in estimating the total uncertainty (E tot ) of the annual carbon balance: The ecosystem WUE was defined as the ratio of GPP to ET, i.e. H 2 O flux: where daily means of GPP and ET were used. The ET data corresponding to a latent heat flux lower than 30 W m −2 were discarded (Abraha et al., 2016). were taken along the 1-m core at 16 points, and soil organic carbon (SOC, kg m −2 ) content in each subsample was analysed using a VarioMax CN analyser (Elementar Analysensysteme GmbH, Germany).

Meteorological conditions
The annual mean air temperature at the study site was 7.6 • C and 7.7 • C in the first and second measurement year, respectively.
Both years were warm compared to the long-term (1981-2010) average of 5.4 • C measured at a nearby weather station (Pirinen et al., 2012). The annual precipitation sum was lower in the first year (473 mm) and higher in the second year (855 mm) than

Fluxes
At the beginning of the measurements, the net CO 2 fluxes were negative (Fig. 4), and the air and soil temperatures were already well above 10 • C (Fig. 5). Net uptake was observed until the first harvest around mid-June 2018. This harvest and the following management events during the that growing season induced large short-term variations in the CO 2 fluxes. Similarly, in the second study year, large impacts on CO 2 fluxes were observed after the management events. During the growing season, the  Seasonal patterns were observed also in the H 2 O fluxes (Fig. 4). In the spring, the ecosystem ET started to increase reaching The experimental field was harvested and fertilised twice during each of the studied growing seasons (Table 1). The effect of management was investigated by comparing the mean fluxes 5 days before and after the harvest dates (Table A1) In the first growing season, the first and second fertilisation events with organic substances increased NEE by 0.27 and 0.08 mg CO 2 m −2 s −1 , respectively, i.e. diminished the CO 2 sink (Fig. 4, Table A1). During the 5 days after the harvest in May 2019, the field acted as a CO 2 source. A similar trend was not observed in June 2019, as mineral fertiliser was used and thus no organic substances were added to the soil. Each of the fertilisation events were followed by rain within the next 265 5 days. However, the mean soil moisture remained either the same or decreased slightly (Fig. 5, Table A1). Furthermore, the mean air temperature increased after the fertilisations in July 2018 and May 2019, potentially affecting CO 2 fluxes. After the fertilisation events with organic substances in July 2018, August 2018 and May 2019, the mean PAR was 7%, 29% and 12% lower, respectively, than the 5-day mean before the fertilisation, complicating the interpretation of fertilisation impacts on the CO 2 fluxes. The effect of management on H 2 O fluxes could not be disentangled from the present data (Fig. 4b). 270 The PI calculated from the flux data was consistent with the seasonal changes in the LAI derived from Sentinel-2 images (Fig. 5d). The higher LAI in 2019 indicated that there was more photosynthesising green biomass before the first and second harvest compared to 2018. The effect of larger leaf area was also observed in the differences in the photosynthetic capacity (GP 1200 ) of the grassland between the study years (Fig. 6a). The years differed significantly (p < 0.05) in terms of GP 1200 at all levels of LAI (>1). Larger LAI values were observed throughout 2019, indicating that grass was growing better than in 2018.

275
Furthermore, the grassland was photosynthesising more efficiently with the same leaf area in 2019 than in the previous year ( Fig. 6a).

Water use efficiency
The ecosystem WUE estimate showed different seasonal variation during the studied growing seasons (Fig. 7). Generally, WUE was higher in 2019 than in 2018 throughout the growing season. WUE increased before the first harvest around mid-280 June in both years, indicating more efficient CO 2 uptake in terms of water use than during the spring. The 5-day mean WUE before the first harvest was 2.6 and 2.9 g CO 2 (kg H 2 O) −1 in 2018 and 2019, respectively. Due to the harvest, it dropped to 0.8 g CO 2 (kg H 2 O) −1 in 2018 and to 2.2 g CO 2 (kg H 2 O) −1 in 2019. During the latter growing season, WUE increased steadily towards 4 g CO 2 (kg H 2 O) −1 until the second harvest in August, whereas WUE remained predominantly below 2 g CO 2 (kg H 2 O) −1 during the same period in 2018. In the end of August and early September, WUE was at the same level in both years.

285
The LAI derived from the Sentinel-2 data was compared to the daily WUE values (Fig. 6b) to further interpret the relation between vegetation status and ecosystem WUE. While WUE was on average lower in 2018 than 2019, the difference at a given LAI was not significant (p > 0.05). However, in both years the daily WUE increased in a similarly linear manner in relation to LAI.

Carbon balance and soil carbon content 290
The carbon balance of the studied grass field was -50 ± 68 g C m −2 yr −1 , i.e. not different from zero, in the first year, while the balance of the second year was negative, -118 ± 24 g C m −2 yr −1 , i.e. the field acted as a net carbon sink (Table 2). All the components of the carbon balance were smaller in the first year than in the second one, GPP by 30%, R eco by 25% and management by 94%.   14 https://doi.org/10.5194/bg-2020-422 Preprint. Discussion started: 30 November 2020 c Author(s) 2020. CC BY 4.0 License. There was a major difference in the CO 2 balances between the growing seasons (Table 3). In 2019, the growing season net 295 uptake was 78%, GPP 49% and R eco 42% higher than in 2018.
The average soil carbon content in the 1-m layer was 16.59 ± 2.25 kg m −2 (average ± standard deviation), with the highest SOC found in the top 30-cm layer (Fig. 8). The carbon balance of 2018 was 0.3% of the average SOC, and in 2019 this ratio was 0.7%. On average, the annual carbon input to the soil accounted for 0.5% of the SOC.  GPP of European managed grasslands to be in the range of -2107 to -1410 g C m −2 yr −1 . Our results fall in to the lower range or below these GPP values. The annual R eco in Qvidja was also varying between the study years (972 and 1291 g C m −2 yr −1 ). Globally, the annual R eco of managed grasslands is reported to vary within a wide range from 31 to 2150 g C m −2 yr −1 . The average R eco was 1445 and 647 g C m −2 yr −1 on the intensively and extensively managed grasslands, respectively (Gilmanov et al., 2010), and the annual R eco in Qvidja falls between these values in both study years. Regarding only the 310 European grasslands, the annual R eco is reported to vary between 494 and 1623 g C m −2 yr −1 (Gilmanov et al., 2007). Our observations are thus also within this range.
The carbon balance of the grass field ecosystem in Qvidja was close to neutral (-50 ± 68 g C m −2 yr −1 ) in the first study year (  . These nine grassland sites acted mainly as net carbon sinks, the annual net carbon balance ranging from -446 to 251 g CO 2 -C eq. m −2 yr −1 , where 13 of the 17 measured annual balances were negative. Our study site falls into this range. In comparison, the Finnish agricultural sites measured so far are generally carbon sources (Heikkinen et al., 2013;Shurpali et al., (two-year average NEE -259 g C m −2 yr −1 ) for a grassland site on mineral soil than we observed in Qvidja. However, by considering the total carbon balance of the system by taking into account the carbon fluxes caused by management, it was concluded that their site acted as a net carbon source. Mineral fertilisers were used during their study, and thereby no carbon was imported to the field to compensate for the biomass removal from the system as harvests. Similar carbon flux patterns related to management were reported by Eichelmann et al. (2016). The annual NEE of the agricultural grassland in their study 325 in Canada was more negative (average NEE -405 g C m −2 yr −1 ) than the NEE in Qvidja. However, the two-year mean annual carbon balance of the Canadian field was positive when biomass removal was taken into account, i.e. the field was a net source of carbon. It is noteworthy that the yield in Qvidja was substantially smaller than at the other two study sites (Lind et al., 2016;Eichelmann et al., 2016), at which the total balance became positive when the management was taken into account. However, as the harvested grass was used as feed for farm animals, there was no need for a higher yield at Qvidja in either of the study 330 years.
Analysis of the weather variables in Qvidja indicated that temperature and moisture conditions were associated with the differences in fluxes and carbon balance between the study years. The growing season was warmer and drier in 2018 than 2019, with 13% lower mean soil moisture, 32% lower precipitation, 2.2 • C higher mean air temperature and 12% higher radiation during the growing season, and substantially smaller fluxes were observed in the first year. This is in accordance 335 with Shurpali et al. (2009) who observed a positive correlation between the uptake of CO 2 and both soil moisture and air temperature on another Finnish agricultural grassland. According to their conclusions, moderate temperature with high soil moisture favoured CO 2 uptake. Furthermore, Flanagan et al. (2002) and Kurc and Small (2007) concluded that photosynthesis of grassland favours rather wet summer conditions. These findings would support the conclusion that low soil moisture and high temperatures were the main factors limiting CO 2 uptake at our study site in the summer 2018.

340
To answer our first research question, we conclude that there were notable year-to-year differences in the carbon balances, but the reason behind this variation remains partly open, as weather, grass age and grass leaf area all showed different dynamics between the study years. In Finland, it is typical to grow grasslands for 3-4 years before the grass renewal. In Qvidja, the grass was not renewed between the study years, which may have led to the larger fluxes observed in the second year when the grass root system, for instance, was possibly more developed enhancing carbon uptake. Furthermore, the leaf area was larger, and of the more developed grass may have increased carbon uptake. The lower leaf area during the first year was most probably also connected to the dry summer, as shortage of water is a growth-limiting factor. Besides the leaf area, the photosynthetic potential per leaf area was lower in the first year, indicating either drought stress or shortage of nutrients, as temperature, a widely limiting factor in northern latitudes, was high enough during both summers not to restrict photosynthesis. In any case, 350 a more specific analysis of the driving and inhibiting environmental factors will require a longer measurement period.
Our second research question concerned the drought-related restrictions of photosynthesis. It has been widely recognised that in dry conditions plants are able to reduce transpiration by stomatal regulation (Pirasteh-Anosheh et al., 2016). However, grasses seem to limit stomatal functions only in severe, prolonged drought conditions (Wolf et al., 2013;Xu et al., 2019), and thus occasional or seasonal drought events may not be observed in the ecosystem WUE of grasslands. In our study, WUE 355 values were predominantly lower in 2018 than in 2019. This was most probably explained by the differences in LAI as the relationship between WUE and LAI was similar in both growing seasons (Fig. 6b). Furthermore, the drier conditions with high temperatures in the summer 2018 may have resulted in a decoupling of assimilation and transpiration and in temperatureinduced downregulation of GPP (Gharun et al., 2020), as ET was similar in both years (Table 3). Therefore, the clearly lower leaf-area-based photosynthetic capacity (GP 1200 ) in 2018 compared to 2019 probably indicates drought related stress in pho-360 tosynthetic processes despite the similar leaf-area-based WUE (Fig. 6). It is noteworthy that the WUE analysis was performed by means of the total ecosystem ET rather than plant transpiration, which would have enabled a more direct determination of the actual plant WUE and thus a simpler interpretation of plant processes and their relation to LAI. In general, WUE at our study site varied between 0 and 5 g C (kg H 2 O) −1 . This is consistent with the WUEs observed for northern grasslands (0-7 g C (kg H 2 O) −1 ) (Tang et al., 2014).

365
The different management practices, such as fertilisation and the choice of grass cutting height, were slightly different in the first and second year, which probably had an impact on the carbon balances. In June 2018, a conventional cutting height of 6 cm was used, whereas in the other harvests the grass was cut at 15 cm. The higher cutting height may have enhanced the regrowth of the grasses, especially in the more favourable weather conditions in 2019, and with a larger leaf area higher CO 2 uptake was observed right after the harvest. Only after the 6-cm harvest, the field turned to a net source of CO 2 . With a low 370 cutting height, it was more likely that the grass was cut below the growing point, particularly in dry conditions, which affects the stand longevity and stress tolerance (Jones and Tracy, 2018). As the weather was warm and dry during the harvest events in June in both years, a higher cutting height may have served as a vital management improvement.
The field was mainly fertilised with organic substances, and thus carbon was imported to the system, affecting the net carbon balance. After each of the fertilisation events with organic material, the respiration of the field seemed to increase, whereas 375 mineral fertilisation was not observed to have an immediate effect on CO 2 fluxes. Increased respiration was likely to occur due to microbial activity of the organic fertilisers. Gilmanov et al. (2007) observed on a Danish agricultural grassland that, although the application of manure increased respiration, also the plant uptake of CO 2 was notably higher than at the other sites studied. Fornara et al. (2016) also concluded, based on their 43-yr study, that manure fertilisation substantially increased soil carbon sequestration of a grassland ecosystem in Northern Ireland. Although the type of the organic fertiliser possibly plays a crucial 380 role, the application of carbon to the system has a direct effect on the carbon balance, but there is also an indirect effect on its components R eco and GPP via soil and plant functions.
Concerning our final research question on the relation of possible carbon sink to the international "4 per 1000" carbon sequestration initiative (Minasny et al., 2017), our results show that, on average, the field acted as a net annual carbon sink by increasing the soil carbon content by 0.5% annually over the studied period. Thus the site fulfilled the goal of the "4 per 1000" 385 initiative and contributed to the short-term climate change mitigation. Furthermore, the annual carbon balance of our second study year (-118 g C m −2 yr −1 ) is in the upper range of annual carbon sequestration potential (80-120 g C m −2 yr −1 ) that is evaluated to be attainable with improved management practices (Lal, 2016). Thus, this study demonstrates the potential for a positive impact of northern agricultural grasslands in terms of climate change mitigation.

Errors and uncertainties
Uncertainties with the data are mainly related to the gaps in the measurement data, which required gap-filling of those periods with modelled data. The length of a gap increases the related uncertainty, but in our data there were only three longer gaps (4, 8 and 9 days), which all occurred during the first winter, when temperatures were low and only minor fluxes could have been observed. All the other gaps were shorter than 3 days. However, each gap contributed to the uncertainty and were included in the carbon balance calculations. Further uncertainties, which were not included in the error estimates, were caused by the the 395 soil temperature modelling for the first study year and the management flux estimates.
Carbon balance was calculated based on the ecosystem-atmosphere CO 2 fluxes and the inputs and outputs of harvest and fertilisation. Thus, no other gaseous carbon compounds, such as methane, were considered. Regina et al. (2007) reported that the annual methane exchange of a Finnish clay soil varied between -0.009 and 0.034 g CH 4 m −2 yr −1 during two years in 2000-2002. Thus, based on this estimate, the possible carbon emission from methane accounts for less than 1% of our annual 400 carbon balance.
Leaching of dissolved carbon and emissions of volatile organic compounds may have had an effect on the annual carbon balance. Leaching of carbon from the agricultural soils is mainly driven by the meteorological and hydrological conditions (Manninen et al., 2018), but it is also affected by soil properties (Don and Schulze, 2008). Large variations in soil moisture and temperature and precipitation may increase the solubility of SOM. Generally, however, clay soils retain carbon better than 405 other soil types. Furthermore, ploughing increases leaching as mineralisation of SOM is enhanced. Depending on precipitation and hydrological and chemical properties of the soil, carbon leaching on grasslands may equal approximately to 25% of the annual carbon balance calculated based on NEE, harvest and fertilisation (Kindler et al., 2011). At our study site, the effect of leaching on carbon balance could be assumed to be fairly small in both summers due to low soil moisture and low precipitation.
On the other hand, during wet periods, the leaching may have had a small effect on carbon balance. However, a more precise 410 carbon balance estimate would require further measurements, including leaching and other carbon-containing gases.

Conclusions
The agricultural grassland site located at Qvidja in southern Finland acted as a net carbon sink during the two years studied.
The carbon balance of the first study year was -50 ± 68 g C m −2 yr −1 and in the second year it was -118 ± 24 g C m −2 yr −1 .
We estimated that on average the grassland exceeded the goal of the "4 per 1000" initiative intending to increase soil carbon 415 content. The data presented here act as a basis for the future studies at this site that focus on the conversion from intensive agricultural practices towards more sustainable agricultural management and its impacts on the GHG fluxes on mineral soils in northern conditions. Further research with longer-term measurements would be needed to evaluate the persistence of carbon sequestration and storage. Longer time series would be also essential to study more closely the causes of the interannual variation of GHG fluxes and carbon and water balances at this site.

Appendix A: Effective phytomass index
The PI was used to refine the gap-filling of GPP, especially in the case of long gaps in the nighttime data, based on which GPP was parameterised. PI reflects the development of LAI but, being derived from the daytime NEE measurements, is more 425 dynamic than LAI and thus describes more precisely the course of the photosynthetic activity of plants (Aurela et al., 2001).
PI was derived from the net ecosystem CO 2 exchange data by selecting fluxes at high PAR levels. The PAR limit was set to 700 µmol m −2 s −1 from March to September and 200 µmol m −2 s −1 from October to February. The assumed respiration, i.e. fluxes when PAR < 20 µmol m −2 s −1 , was subtracted from from the NEE data. This was followed by averaging NEE and R eco within a moving window, which was set to 3 days and increased to 5 or 7 days if necessary. Averaging was limited to the 430 harvest dates by decreasing the window size step-by-step to 1.5 days, and similarly increasing it after the harvest. An average GPP was then calculated by subtracting R eco from NEE and normalising to unity to obtain PI. Daily PI values were used for calculating the GPP fluxes. Due to the scarcity of respiration data in July in both years and of the daytime data in winter, linear interpolation was applied to cover the missing daily PI values.
Appendix B: Gap-filling of CO 2 fluxes

435
The flux data set was separated into sections at the harvest dates, and gap-filling was done separately for these sections by first calculating R eco and then GPP. The parameter R 0 was determined for each day from the nighttime data (PAR < 20 µmol m −2 s −1 ) with a 7-day moving window. If there were less than 24 measurements within the time window, its length was increased by 1 day both in the beginning and at the end until enough data were obtained. R 0 was allowed to vary between 0.001 and 1 mg m −2 s −1 . Similarly, the same minimum number of observations and a 3-day moving window was used for determining α 440 and GP max from the observed NEE from which the estimated R eco had been subtracted. α and GP max were allowed to vary between -0.1 and 0 mg µmol −1 , and -5.0 and 0 mg m −2 s −1 , respectively. From 5 December 2018 to 26 March 2019 and from 26 November 2019 to 15 March 2020, with no significant CO 2 uptake, a 5-day moving average was used to fill the gaps in NEE.
Appendix C: Gap-filling of energy fluxes 445 The gaps in the net radiation (R n ) time series were filled with the monthly mean diurnal cycles. Soil heat flux (G) was not measured at our site, so it was estimated from the energy balance closure during the periods when the other energy fluxes were known. Gap-filling of G was done by assuming a constant ratio between G and R n (Liebethal and Foken, 2007). The ratio of 0.24 was calculated with linear regression from the daytime data (between 10:00-15:00). The sensible and latent heat fluxes (Q H and Q E , respectively) were gap-filled based on the procedure described by Kowalski et al. (2003). The gaps in the Table A1. Mean flux and meteorological conditions 5 days before and after management. The management day is not included.