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

A significant proportion of the global carbon emissions to the atmosphere originate from agriculture. Therefore, continuous long-term monitoring of CO₂ fluxes is essential to understand the carbon dynamics and balances of different agricultural sites. Here we present results from a new eddy covariance flux measurement site located in southern Finland. We measured CO₂ and H₂O fluxes at this agricultural grassland site for 2 years, from May 2018 to May 2020. In particular the first summer experienced prolonged dry periods, which affected the CO₂ fluxes, and substantially larger fluxes were observed in the second summer. During the dry summer, leaf area index (LAI) was notably lower than in the second summer. Water use efficiency 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₂ fluxes and management measures, which included input of carbon as organic fertilizers and output as yield. The carbon balance of the field was −57 ± 10 and −86 ± 12 g C m⁻² yr⁻¹ in the first and second study years, respectively.

1 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., 2009, 2017; Lal, 2016; Paustian et al., 2000; Smith, 2008). Currently, agriculture is responsible for more than 10 % of the global 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₂) emissions from land use (Ceschia et al., 2010; Reinsch et al., 2018; Yang et al., 2019). 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 biodiversity and cover cropping, would not only help to mitigate climate change but also to restore soil quality and fertility. In particular managed grasslands as part of agricultural systems have a high potential for substantial soil carbon sequestration (Soussana et al., 2010; Gilmanov et al., 2010; Yang et al., 2019). The importance of increasing soil organic carbon (SOC) content of agricultural soils has recently attained more attention, and the “4 per mille Soils for Food Security and Climate” initiative was launched at the 21st Conference of the Parties to the United Nations Framework Convention on Climate Change in Paris in 2015 (Minasny et al., 2017). The aim of this initiative is to increase the soil carbon stock on all land surfaces in the upper 2 m on average by 0.4 % annually. The possible increase in carbon content is largely dependent on the soil properties, e.g. clay content (Johannes et al., 2017; Minasny et al., 2017). This would be enough to sequester carbon from the atmosphere by an amount equivalent to the annual anthropogenic GHG emissions. However, the initiative states that the most potential SOC increases can be achieved on managed agricultural lands. In that case, the “4 per 1000” means increasing SOC at the top 1 m layer of agricultural soils by 0.4 % annually. That would effectively offset approximately 20 %–35 % of the global GHG emissions.

Agricultural ecosystems are highly prone to impacts of climate change, which induces a risk for food production. One of the possible impacts of climate change on agricultural ecosystems is associated with the changes in seasonal weather conditions and the resulting alteration in the carbon and water balance of these ecosystems (Ciais et al., 2014; Donnelly et al., 2017; Harrison et al., 2019). Severe drought events and storms causing considerable damage to agriculture have already been observed across Europe (Ciais et al., 2005; Wolf et al., 2013; Bastos et al., 2020). Moreover, adverse climatic impacts may be amplified by current and prior land use practices if they have not supported ecosystem resilience (Brunsell et al., 2014). For instance, a deeper root system is likely to buffer the negative impacts of climate variability. Also, high plant species diversity, compared to monocultures, favours the efficiency of plant water consumption and resilience to drought (De Boeck et al., 2006). As gross primary production (GPP) is closely related to ecosystem evapotranspiration (ET) via stomatal functions (Fricker and Willmer, 2012), changes in terrestrial water balance are potentially reflected in GPP and thus in the carbon balance of agricultural grasslands. The effect of water stress can be studied, for instance, by analysing ecosystem water use efficiency (WUE), i.e. the amount of carbon assimilated per unit of water lost by transpiration (Steduto, 1996). Generally, the productivity of a grassland ecosystem correlates with WUE, and thus ecosystems with a high productivity usually also have a high WUE (Hu et al., 2008). Environmental factors regulate WUE via effects on stomatal conductance and GPP, and during prolonged drought periods, for example, temperature-induced downregulation of GPP may reduce WUE of grasslands in particular (Gharun et al., 2020). Furthermore, the WUE response depends on the intensity of the drought (Xu et al., 2019). However, the drought effects are also strongly related to season, as Wolf et al. (2013) reported that the WUE of Swiss grassland ecosystems did not respond to a spring drought and Bastos et al. (2020) concluded that the spring weather may either amplify or dampen the carbon and water dynamics during the following summer.

Better understanding of climatic impacts of agriculture and the effects of improved practices from the perspective of soil health and vitality is needed in order to develop tools for improved environmental management of these ecosystems. Continuous long-term measurements of the atmosphere–ecosystem fluxes are needed to identify the key factors affecting carbon dynamics of different ecosystems, to quantify the resulting carbon balance and its components, and to verify soil carbon and ecosystem models. Moreover, high-quality GHG flux data are needed for a reliable, global measuring, reporting and verification system of agricultural carbon fluxes and soil carbon sequestration and stability (Smith et al., 2020).

The eddy covariance (EC) method is widely used for measuring CO₂ and energy fluxes in different ecosystems and climatic conditions (Aubinet et al., 2012). The high-frequency measurements provided by EC allow a direct quantification and analysis 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 northern European countries (Shurpali et al., 2009; Lind et al., 2020; Jensen et al., 2017).

The aim of this study is to investigate, based on EC measurements, CO₂ exchange between the atmosphere and a managed forage grassland in southern Finland. In particular, we had three specific research questions:

1. What is the magnitude of the annual carbon balance and its components?

2. Does the grass photosynthesis indicate occasional drought-related responses?

3. How does the possible carbon sink relate to the carbon sequestration objective of the “4 per 1000” initiative?

For the purposes of this study, we collected field data on the net exchange of CO₂ and H₂O, soil and vegetation properties, and meteorological variables on an agricultural grassland in southern Finland for 2 years, from May 2018 to May 2020.

2 Material and methods

2.1 Site description

The flux measurements were conducted at the Qvidja farm in southern Finland (60.29550^∘ N, 22.39281^∘ E; elevation 5 m) from 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 more sustainable and environmentally friendly farming by increasing the use of organic fertilizers and perennials, restricting the use of pesticides and increasing plant species biodiversity. The current grass and clover mixture was sown as an undergrown species with broad bean in spring 2017. The predominant species were timothy (Phleum pratense), meadow fescue (Festuca pratensis) and white clover (Trifolium repens).

Figure 1 Experimental field with the sectors representing the target area that covers 3.9 ha. Eddy covariance tower is located in the centre of the sectors. EC data from wind directions from 30 to 140^∘ were discarded due to another experimental plot located in that part of the field. (Orthophoto from National Land Survey of Finland.)

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 to drought, oversowing 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 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 on 23 September. In 2019, the grass was harvested on 11 June and 20 August. In June 2018, a conventional cutting height of 6 cm was used, whereas in the other harvests the grass was cut at 15 cm.

In 2018, the field was fertilized twice, on 16 July and 24 August, with 2800 and 1800 kg ha⁻¹ of NK-molasses, respectively (Table 1). NK-molasses is 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⁻¹ of organic carbon. In addition, it contained potassium and small proportions of sulfur, magnesium, calcium and sodium.

Table 1 Different management events and their C inputs (fertilization) and C outputs (harvest). During the cutting in August 2018, the grass was not collected and thus did not result in any C flux allocated to management.

In May 2019, the field was fertilized with a mixture of side products from industries of starch potato processing, biowaste processing and ethanol production out of sawdust. This fertilization mixture contained 70 % (of dry weight) of OM, 1.3 % of nitrogen, 0.2 % of phosphorus, 3 % of potassium and 0.4 % of sulfur, as well as small amounts of calcium, magnesium, zinc, copper and manganese. Approximately 4600 kg ha⁻¹ was applied to the field on 8 May (Table 1). On 26 June after the first harvest, 220 kg ha⁻¹ of mineral fertilizers was applied. This fertilizer contained 23 % of nitrogen, 10 % of phosphorus and 8 % of potassium.

2.2 Measurement setup

The CO₂ and H₂O fluxes were measured with the micrometeorological EC method. The flux measurements started on 3 May 2018, and here we analysed data collected from 4 May 2018 to 3 May 2020. From this point on, the periods of 4 May 2018–3 May 2019 and 4 May 2019–3 May 2020 are referred to as the first and second EC measurement years, respectively.

The EC instrumentation consisted of an enclosed infrared gas analyser (LI-7200, LI-COR Biosciences, NE, USA), which detects the CO₂ and H₂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⁻¹, and the length of the 4 mm stainless-steel inlet tube with 2 µm Swagelok sinter was 0.8 m. The CO₂ measurements were regularly checked with zero and span gases, and the LI-7200 was recalibrated when necessary. The H₂O measurements were compared with the data obtained from a dedicated humidity sensor; no recalibration was necessary.

The micrometeorological sign convention is used throughout the paper, with a negative value indicating 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 profiles at the depths of 5, 10 and 30 cm (Pt100 IKES sensors, Nokeval Oy, Nokia, Finland). The soil temperature data were collected with a Vaisala QML201C datalogger (Vaisala Oyj, Vantaa, Finland). Photosynthetically active radiation (PQS PAR sensor, Kipp & Zonen B.V., Delft, the Netherlands), global and reflected solar radiation (CMP3 radiometer, Kipp & Zonen), and air temperature and relative humidity (Humicap HMP155, Vaisala Oyj) were measured at the height of 1.8 m. In addition, precipitation was measured with a weighing rain gauge (Pluvio2, OTT HydroMet GmbH, Kempten, Germany). Meteorological measurements started on 8 May 2018, and the data were recorded as 30 min averages, excluding the precipitation which was recorded as 1 min values. Snow depth was recorded at the weather station of Kaarina Yltöinen.

The leaf area index (LAI) data were obtained from the Sentinel-2 satellite as daily values on the clear-sky days. LAI was calculated from the Sentinel-2 bottom-of-atmosphere products (L2A) using the Google Earth Engine (GEE) and a Python implementation of the Biophysical Processor toolbox (Weiss and Baret, 2016) available in Sentinel Application Platform (SNAP) software. The cloudy, cloud-shadowed and snowy data were filtered out using the scene classification band available in the L2A products.

2.3 Eddy covariance data processing

The turbulent fluxes were determined as the covariance between the variations in 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 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₂ mixing ratio was set to 410 ppm in these calculations.

The following acceptance criteria were applied to screen the 30 min averaged CO₂ flux data: number of spikes in the raw data <150 of 18 000, relative stationarity of CO₂ flux (Foken et al., 2012) <50 %, mean CO₂ mixing ratio >380 ppm, variance of CO₂ mixing ratio <15 ppm² between April and September and <5 ppm² 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 Meixner (2001). For these, we applied a friction velocity limit of 0.06 m s⁻¹ and a cumulative footprint limit of 0.7. The further screening applied to H₂O fluxes included H₂O flux >0, relative stationarity of H₂O flux <50 % and variance of H₂O mixing ratio <1 (mmol mol⁻¹)². After applying these filtering criteria, the coverage of CO₂ and H₂O flux data accepted for further analysis was 44 % and 30 % of all the 30 min periods (i.e. total of 35 088 time steps) during the 2 measurement years, respectively (for CO₂, day/night 55 %/33 %, April–September/October–March 48 %/38 %; for H₂O, day/night 49 %/11 %, April–September/October–March 41 %/16 %). Most of the accepted CO₂ and H₂O flux data were collected when the wind direction was in the south-southwest sector (Fig. 2).

Figure 2 Number of accepted flux measurements within 20^∘ sectors around the flux tower during the first and second years. Data from 30 to 140^∘ were discarded.

2.4 Flux partitioning and gap-filling

To calculate CO₂ balances and to analyse the components of the net ecosystem exchange between the field and the atmosphere, the measured CO₂ 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:

$$\begin{array}{}\text{(1)}& \mathrm{NEE}=\mathrm{GPP}+R_{\mathrm{eco}}.\end{array}$$

The gap-filled GPP and R_eco were calculated with empirical response functions by first fitting these functions to the flux data. R_eco was expressed as a function of temperature (Lloyd and Taylor, 1994):

$$\begin{array}{}\text{(2)}& R_{\mathrm{eco}}=R_{\mathrm{0}}{e}^{E_{\mathrm{0}}\left(\frac{\mathrm{1}}{T_{\mathrm{1}}}-\frac{\mathrm{1}}{T_{\mathrm{a}}-T_{\mathrm{0}}}\right)},\end{array}$$

where R₀ is the respiration rate (mg m⁻² s⁻¹) at the reference temperature of 283.15 K, T₀=227.13 K, T₁=56.02 K, E₀ is ecosystem sensitivity coefficient (Lloyd and Taylor, 1994) that describes the temperature response of soil respiration and T_a is the air temperature.

GPP was modelled as a function of photosynthetically active radiation (PAR, µmol m⁻² s⁻¹) as

$$\begin{array}{}\text{(3)}& \mathrm{GPP}={\displaystyle \frac{\mathit{\alpha }\times \mathrm{PAR}\times {\mathrm{GP}}_{\max }}{\mathit{\alpha }\times \mathrm{PAR}+{\mathrm{GP}}_{\max }}},\end{array}$$

where α is the apparent quantum yield (mg µmol⁻¹), and GP_max denotes the asymptotic CO₂ uptake rate in optimal light conditions (mg m⁻² s⁻¹). Further details on the gap-filling procedure are provided in Appendix A. Energy fluxes were gap-filled following the description in Appendix B.

To study the differences in photosynthetic capacity of the grass field between the two growing seasons, daily GP₁₂₀₀ values were calculated with the estimated α and GP_max values; i.e. GPP was normalized to PAR = 1200 µmol m⁻² s⁻¹.

2.5 Net ecosystem carbon balance

In this study, the system boundaries include the main components of the carbon balance of the field ecosystem studied. The carbon balance was calculated by adding up the 30 min NEE fluxes, the imported carbon in the form of organic fertilizers and the carbon removed as harvested biomass:

$$\begin{array}{}\text{(4)}& C_{\text{balance}}=C_{\mathrm{H}}+C_{\mathrm{F}}+\sum _{i=\mathrm{1}}^n{\mathrm{NEE}}_i,\end{array}$$

where C_H is the amount of carbon in harvested biomass, C_F is the amount of carbon in imported fertilization and n is the total number of time steps in the period for which the balance was calculated. Thus, the carbon balance indicates the net ecosystem carbon balance as defined by Chapin et al. (2006) without the contribution of carbon monoxide, methane, volatile organic and particulate compounds, or leaching. This balance is commonly called the net biome production (Kutsch et al., 2010). 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.

2.6 Uncertainty analysis

The CO₂ balance between the field and the atmosphere, which is calculated based on the EC measurements, includes multiple potential error sources. Uncertainties 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 (Aubinet et al., 2012). 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), were included in the uncertainty estimate:

$$\begin{array}{}\text{(5)}& E_{\mathrm{meas}}=\sqrt{\sum _{i=\mathrm{1}}^n({\mathrm{NEE}}_{\mathrm{meas},i}-{\mathrm{NEE}}_{\mathrm{mod},i}{)}^{\mathrm{2}}},\end{array}$$

where NEE_meas is the filtered 30 min flux, NEE_mod is the corresponding modelled NEE (Eqs. 1–3) and n is the number of measured data.

$$\begin{array}{}\text{(6)}& E_{\mathrm{gap}}=\sqrt{\sum _{i=\mathrm{1}}^N(E_{\mathrm{GPP},i}^{\mathrm{2}}+E_{R_{\mathrm{eco}},i}^{\mathrm{2}})},\end{array}$$

where E_GPP and $E_{R_{\mathrm{eco}}}$ are the errors of modelled GPP and R_eco, respectively. N is the number of gaps in the data.

The standard error propagation principle was used in estimating the total uncertainty (E_tot) of the annual carbon balance:

$$\begin{array}{}\text{(7)}& E_{\mathrm{tot}}=\sqrt{E_{\mathrm{meas}}^{\mathrm{2}}+E_{\mathrm{gap}}^{\mathrm{2}}}.\end{array}$$

2.7 Water use efficiency

The ecosystem WUE was defined as the ratio of GPP to ET, i.e. H₂O flux:

$$\begin{array}{}\text{(8)}& \mathrm{WUE}={\displaystyle \frac{\mathrm{GPP}}{\mathrm{ET}}},\end{array}$$

where daily means of GPP and ET were used. The ET data corresponding to a latent heat flux lower than 30 W m⁻² were discarded (Abraha et al., 2016). Furthermore, days with precipitation were eliminated in order to obtain a signal that is dominated by transpiration.

2.8 Soil carbon storage

Soil carbon content was determined from 1 m deep core samples taken within the flux source area. The samples were taken in October 2018 using a hydraulic corer installed on a tractor. The diameter of the sample cylinder was 151 mm. Subsamples were taken along the 1 m core at 16 points, and soil organic carbon (SOC, kg m⁻²) content in each subsample was analysed using a VarioMax CN analyser (Elementar Analysensysteme GmbH, Germany).

3 Results

3.1 Meteorological conditions

The annual mean air temperature at the study site was 7.6 and 7.7 ^∘C in the first and second measurement years, 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 the long-time average (679 mm).

The thermal growing season, defined here as the period when the daily mean temperature permanently exceeded 5 ^∘C, started on 14 April in 2018, i.e. before the EC measurements started. In 2019 and 2020, the thermal growing season began on 16 and 18 April, respectively. The thermal growing season ended on 17 November and 26 October in 2018 and 2019, respectively. Thus, the thermal growing season length was 218 d in 2018 and 194 d in 2019. Meteorological conditions during the main growing season between May and September varied substantially between the 2 years. The mean air temperature during these periods was 16.7 and 14.5 ^∘C in 2018 and 2019, respectively. During the same period, the mean daily PAR was about 12 % higher in 2018 than in 2019 (460 vs. 410 µmol m⁻² s⁻¹), while the precipitation sum was 32 % lower (212 vs. 312 mm).

During winter 2018–2019, permanent snow cover was recorded from 17 December 2018 to 26 March 2019. In the following winter (2019–2020), there were only two short snow-cover periods: 5–8 February and 30–31 March 2020. The maximum snow depth in the first winter was 33 cm, whereas in the second winter it was 3 cm. The mean wintertime (November–March) air temperature was −0.2 ^∘C in 2018–2019 and 2.2 ^∘C in 2019–2020.

Soil moisture content at the depth of 10 cm varied between 0.16 and 0.55 m³ m⁻³ during the study period. On several occasions, the daily mean soil moisture dropped to about 0.2 m³ m⁻³. During the growing seasons, such low values indicate substantial drought, while in the winter, rapid data drops were likely related to soil freezing. The average soil moisture during the growing season in 2019 was higher than in 2018 (0.30 vs. 0.26 m³ m⁻³). As a result of the higher precipitation in 2019, soil moisture occasionally increased up to 0.4 m³ m⁻³, i.e. close to the saturated values observed in winter.

3.2 CO₂ and H₂O fluxes

At the beginning of the measurements, the net CO₂ fluxes were negative (Fig. 3), and the air temperature was already well above 10 ^∘C (Fig. 4). Net uptake was observed until the first harvest around mid-June 2018. This harvest and the following management events during that growing season induced large short-term variations in the CO₂ fluxes. Similarly, in the second study year, large impacts on CO₂ fluxes were observed after the management events. During the growing season, the mean NEE was −0.13 and −0.21 mg CO₂ m⁻² s⁻¹ in 2018 and 2019, respectively. During the wintertime, no significant CO₂ uptake occurred, and the positive fluxes were small compared to the nocturnal fluxes in summer. The mean measured NEE between December 2018 and February 2019 was 0.03 mg CO₂ m⁻² s⁻¹, and during the same period in 2019–2020 it was 0.04 mg CO₂ m⁻² s⁻¹.

Figure 3 Accepted 30 min (a) net ecosystem exchange (NEE) and (b) H₂O flux measurements from May 2018 to May 2020. Vertical lines with H and F indicate harvest and fertilization, respectively.

Seasonal patterns were also observed in the H₂O fluxes (Fig. 3). In the spring, the ecosystem ET started to increase, reaching the highest levels between June and August, after which it gradually decreased to wintertime values, i.e. close to zero. The mean growing season H₂O flux was 34.7 mg H₂O m⁻² s⁻¹ in 2018 and 35.5 mg H₂O m⁻² s⁻¹ in 2019. The wintertime (December–February) mean H₂O flux was 3.6 and 3.7 mg H₂O m⁻² in 2018–2019 and 2019–2020, respectively.

The experimental field was harvested and fertilized twice during each of the studied growing seasons (Table 1). The effect of management was investigated by comparing the mean fluxes 5 d before and after the harvest dates (Table C1). The harvest in June 2018 changed the mean CO₂ flux from a net sink of −0.28 mg CO₂ m⁻² s⁻¹ to a source of 0.03 mg CO₂ m⁻² s⁻¹, i.e. increased the net efflux by 0.31 mg CO₂ m⁻² s⁻¹. The first harvest of 2019 increased NEE by 0.47 mg CO₂ m⁻² s⁻¹, but as the pre-harvest mean NEE was −0.50 mg CO₂ m⁻² s⁻¹, the field remained as a net sink. As a result of the second harvest on 23 September 2018, the mean sink decreased from −0.10 to −0.02 mg CO₂ m⁻² s⁻¹, while the harvest on 20 August 2019 caused the sink to change from −0.25 to −0.02 mg CO₂ m⁻² s⁻¹. Thus, after all the harvests with a cutting height of 15 cm, the mean sink rate was diminished to −0.02 or −0.03 mg CO₂ m⁻² s⁻¹.

In the first growing season, the first and second fertilization events with organic substances increased NEE by 0.27 and 0.08 mg CO₂ m⁻² s⁻¹, respectively, i.e. diminished the CO₂ sink (Fig. 3, Table C1). During the 5 d after the harvest in May 2019, the field acted as a CO₂ source. A similar trend was not observed in June 2019, as mineral fertilizer was used, and thus no organic substances were added to the soil. Each of the fertilization events were followed by rain within the next 5 d. However, the mean soil moisture at the depth of 10 cm either remained the same or decreased slightly (Fig. 4, Table C1). Furthermore, the mean air temperature increased after the fertilizations in July 2018 and May 2019, potentially affecting CO₂ fluxes. After the fertilization 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 d mean before the fertilization, complicating the interpretation of fertilization impacts on the CO₂ fluxes. The effect of management on H₂O fluxes could not be disentangled from the present data (Fig. 3b).

Figure 4 Daily mean (a) air and soil (depth = 0.05 m) temperature, (b) photosynthetically active radiation (PAR), (c) precipitation and soil moisture (depth = 0.1 m), (d) leaf area index (LAI), and (e) daily mean NEE, GPP, R_eco and cumulative carbon flux from May 2018 to May 2020.

The LAI derived from Sentinel-2 images (Fig. 4d) varied greatly between the years. The higher LAI in 2019 indicated that there was more photosynthesizing green biomass before the first and second harvests compared to 2018. The effect of larger leaf area was also observed in the differences in the photosynthetic capacity (GP₁₂₀₀) of the grassland between the study years (Fig. 5a). The years differed significantly (p<0.05) in terms of GP₁₂₀₀ at all levels of LAI (>1). Larger LAI values were observed throughout 2019, indicating that grass was growing better than in 2018. Furthermore, the grassland was photosynthesizing more efficiently with the same leaf area in 2019 than in the previous year (Fig. 5a).

Figure 5 (a) Daily photosynthetic capacity (GP₁₂₀₀) and (b) water use efficiency (WUE) as a function of leaf area index (LAI) during the two growing seasons. Grey areas represent the uncertainty bands.

3.3 Water use efficiency

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

The LAI derived from the Sentinel-2 data was compared to the daily WUE values (Fig. 5b) to further cast light on the relationship between vegetation status and 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.

Figure 6 Daily water use efficiency (WUE) during two growing seasons.

3.4 Carbon balance and soil carbon storage

The carbon balance of the studied grass field was −57 ± 10 g C m⁻² yr⁻¹ in the first year, and the balance of the second year was −86 ± 12 g C m⁻² yr⁻¹, i.e. the field acted as a net carbon sink in both years (Table 2). The magnitude of all components of the carbon balance were smaller in the first year than in the second one, GPP by 29 %, R_eco by 23 % and management by 96 %. The components in the mean annual CO₂ fluxes between the field and the atmosphere also indicated major differences between the growing seasons (Table 3). In 2019, the magnitude of the growing season NEE was 78 %, GPP 49 % and R_eco 42 % higher than in 2018.

Table 2 The annual carbon balances and their components (g C m⁻² yr⁻¹) for the 2 measurement years. Negative values indicate C input into the ecosystem, whereas positive values indicate C loss. Management (M) is the sum of the C fluxes due to harvest (positive) and fertilization (negative) events (Table 1). The values after ± represent the uncertainty in NEE.

Table 3 Net ecosystem exchange of CO₂ (NEE, g CO₂ m⁻²), its components gross primary production (GPP) and ecosystem respiration (R_eco), and evapotranspiration (ET, mm) during the growing season (4 May to 30 September) in 2018 and 2019.

The average soil carbon storage in the 1 m layer was 16.59 ± 2.25 kg m⁻² (average ± standard deviation), with the highest SOC found in the top 30 cm layer (Fig. 7). To estimate the increase in soil carbon storage, it was assumed that the magnitude of the net carbon balance represented the amount of carbon accumulated in the soil. Furthermore, to evaluate whether the field had the potential to fulfil the “4 per 1000” initiative, the annual net carbon balance was compared to the average soil carbon content. Thus, the estimated increase in soil carbon storage was 0.3 % and 0.5 % in the first and second years, respectively. On average, the annual carbon input to the soil accounted for 0.4 % of the SOC.

Figure 7 (a) Soil organic carbon (SOC) content at different depths in the 1 m deep soil samples, and (b) the total SOC in the samples. Numbers from 1 to 5 indicate sample numbers.

4 Discussion

4.1 Fluxes and carbon balance

There is an urgent need to find evidence-based climate-friendly practices in agriculture in the boreal region, where the growing season is short and varieties differ from those cultivated in the temperate region. The carbon fluxes we measured on the agricultural grassland at the Qvidja farm in southern Finland clearly indicated that this site was a sink of atmospheric carbon. The annual NEE was −62 g C m⁻² yr⁻¹ in the first study year (4 May 2018–3 May 2019) and −216 g C m⁻² yr⁻¹ in the second year (4 May 2019–3 May 2020). The GPP showed notable variation between the study years as the annual GPP was −1121 and −1583 g C m⁻² yr⁻¹ in the first and second years, respectively. Gilmanov et al. (2010) reported a range of −2107 to −1410 g C m⁻² yr⁻¹ for the GPP of European managed grasslands. Our results fall below or in the lower end of this range. The annual R_eco in Qvidja also varied between the study years (1053 and 1362 g C m⁻² yr⁻¹). The annual R_eco in the European grasslands is reported to vary between 494 and 1623 g C m⁻² yr⁻¹ (Gilmanov et al., 2007). Our observations are thus also within this range.

To answer our first research question, we concluded that the carbon balance was negative in both study years (−57 ± 10 and −86 ± 12 g C m⁻² yr⁻¹), and thus the field acted as a net carbon sink during the study period. In comparison, the Finnish agricultural fields measured so far were generally carbon sources when ecosystem–atmosphere CO₂ fluxes, harvests and the carbon supplied to the system as fertilizers were considered (Heikkinen et al., 2013; Shurpali et al., 2009; Lind et al., 2016; Lohila et al., 2004). Lind et al. (2016) reported a slightly more negative annual NEE (2-year average NEE −259 g C m⁻² yr⁻¹) 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 related to biomass removal as grass yield, it was concluded that their site acted as a net carbon source. Mineral fertilizers 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 management-related carbon flux patterns were observed by Eichelmann et al. (2016), who reported a more negative NEE (average −405 g C m⁻² yr⁻¹) for an agricultural grassland in Canada than the NEE in Qvidja; however, the 2-year mean annual carbon balance was positive when biomass removal was taken into account, i.e. the Canadian 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 activities, i.e. harvests and fertilization, were taken into account. The total carbon balance of the field depends greatly both on the amount of organic matter imported to the system as fertilizers and on the harvest yields, which are affected, for instance, by the applied cutting height.

Analysis of the weather variables in Qvidja indicated that temperature and moisture conditions were associated with the differences in CO₂ flux dynamics 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 mean radiation during the growing season, and substantially smaller fluxes were observed in the first year. This is in accordance with Shurpali et al. (2009), who observed a positive correlation between the uptake of atmospheric CO₂ (GPP) and both soil moisture and air temperature on another Finnish agricultural grassland. According to their conclusions, moderate temperature with high soil moisture enhanced CO₂ uptake. Furthermore, Flanagan et al. (2002) and Kurc and Small (2007) concluded that rather wet summer conditions favoured photosynthetic activity in grasslands. These findings would support the conclusion that low soil moisture and high temperatures were the main factors limiting CO₂ uptake at our study site in the summer 2018. However, this question remains partly open, as weather conditions, 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 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 likely to be more developed, enhancing water and nutrient availability and thus reducing the effect of drought stress. Furthermore, the leaf area was larger, and other capabilities, such as microbial symbioses (e.g. de Vries et al., 2020; Harman and Uphoff, 2019; Moreau et al., 2019), of the more developed grass may have increased carbon uptake. The lower leaf area during the first year was most probably also due 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, 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 recognized that in dry conditions plants are able to reduce transpiration by stomatal regulation (Willmer and Fricker, 1996). 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 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 during both growing seasons (Fig. 5b). Furthermore, the drier conditions with high temperatures in the summer 2018 may have resulted in a decoupling of assimilation and transpiration and in temperature-induced 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₁₂₀₀) in 2018 compared to 2019 probably indicates drought-related stress on photosynthetic processes despite the similar leaf-area-based WUE (Fig. 5). 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. Nevertheless, days with even slight precipitation were eliminated from the analysis, and therefore we can assume that during the growing season most of the water flux arises from transpiration. In general, WUE at our study site varied mainly between 0 and 4 g C (kg H₂O)⁻¹. This is consistent with the WUEs observed on northern grasslands (0–7 g C (kg H₂O)⁻¹) (Tang et al., 2014).

The different management practices, such as fertilization and the choice of grass cutting height, were slightly different in the first and second years, 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 grass, especially in the more favourable weather conditions in 2019, and with a larger leaf area higher CO₂ uptake was observed right after the harvest. Only after the 6 cm harvest did the field turn to a net source of CO₂. With a low 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 fertilized with organic substances, and thus carbon was imported to the system, affecting the net carbon balance. After each of the fertilization events with organic material, the respiration of the field increased, whereas mineral fertilization was not observed to have an immediate effect on CO₂ fluxes. Increased respiration was likely to occur due to microbial activity of the organic fertilizers. Gilmanov et al. (2007) observed on a Danish agricultural grassland that, although the application of manure increased respiration, the plant uptake of CO₂ was notably higher than at the other sites studied. Fornara et al. (2016) also concluded, based on their 43-year study, that manure fertilization substantially increased soil carbon sequestration of a grassland ecosystem in Northern Ireland. Although the type of the organic fertilizer possibly plays a crucial 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 the carbon balance to the international “4 per 1000” carbon sequestration initiative (Minasny et al., 2017), our results show that the Qvidja field acted as an annual net carbon sink and had the potential to fulfil the goal of this initiative and to contribute to the short-term climate change mitigation. By considering the carbon balance by accounting for the ecosystem–atmosphere CO₂ fluxes and the carbon fluxes caused by management activities and comparing that to the measured soil carbon content, the carbon storage of the field increased on average by 0.4 % annually over the studied period. To draw a more reliable conclusion about the carbon sequestration, leaching and other carbon-containing compounds must also be considered in further studies about the carbon balance. Furthermore, the number of soil carbon samples should be increased for a more accurate evaluation of the soil carbon storage of the field, even though the variation among the present samples was small. However, the estimated annual carbon balance of our second study year (−86 g C m⁻² yr⁻¹) with improved management practices was within the range of annual carbon sequestration potential (80–120 g C m⁻² yr⁻¹) 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.

4.2 Errors and uncertainties

Uncertainties in the results are mainly related to the gaps in the measurement data, which required gap-filling of missing measurements 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 d), 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 d. 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 involved in the yield measurements and fertilization input estimates, as well as in the fairly scarce sample size of the soil carbon content measurements.

Carbon balance was calculated based on the ecosystem–atmosphere CO₂ fluxes and the inputs and outputs of harvest and fertilization. Thus, no other gaseous carbon compounds, such as methane, were considered. Regina et al. (2007) reported that the annual methane balances of a Finnish clay soil for 2 years were −0.009 and 0.034 g CH₄ m⁻² yr⁻¹. Based on this estimate, the possible carbon emission as methane accounts for less than 1 % of our annual 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 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 other soil types. Furthermore, ploughing increases leaching as mineralization of SOM is enhanced. Depending on precipitation and hydrological and chemical properties of the soil, carbon leaching on grasslands may equal approximately 25 % of the annual carbon balance calculated based on NEE, harvest and fertilization (Kindler et al., 2011). At our study site, the effect of leaching on the annual carbon balance could be assumed to be fairly small in both summers due to low soil moisture and low precipitation. In winter, the leaching may have caused a temporary contribution to the carbon balance during wet periods and thus reduced the increase in soil carbon storage. For a more accurate carbon balance estimate of this site, however, the contribution of leaching and all carbon-containing gases (e.g. CH₄) should be measured and the number of soil carbon storage measurements increased.

5 Conclusions

The agricultural grassland site located at Qvidja in southern Finland acted as a net carbon sink during the 2 years studied. The carbon balance of the first study year was −57 ± 10 g C m⁻² yr⁻¹, and in the second year it was −86 ± 12 g C m⁻² yr⁻¹. When CO₂ fluxes and carbon fluxes caused by management activities were solely accounted for, the soil carbon storage was assumed to have increased by 0.3 % and 0.5 % in 2018 and 2019, respectively, indicating that northern agricultural grasslands have a potential to contribute to climate change mitigation. The data and results presented here act as a basis for the future studies that focus on the conversion of this farm from intensive agricultural practices towards more sustainable agricultural management, especially on the impacts of such a conversion on the GHG fluxes occurring on mineral soils in northern conditions. Even though we could quantify the sink capacity of the field, further research with longer-term measurements is needed to evaluate the persistence of carbon sequestration and storage, and wider measurements of carbon balance components were to be included. Longer time series and broader GHG flux measurements are also essential to more closely study the causes of the interannual variation in GHG fluxes and carbon and water balances at this site, for which the present study provides a baseline.

Appendix A: Gap-filling of CO₂ fluxes

The flux data set was separated into sections at the harvest dates, and gap-filling was done separately for these sections by first parameterizing and calculating R_eco and then GPP. The parameter R₀ was determined for each day from the nighttime data (PAR <20 µmol m⁻² s⁻¹) with a 7 d moving window. E₀ was determined within the same moving window as R₀. If there were fewer than 24 measurements within the time window, its length was increased by 1 d at both the beginning and end until enough data were obtained. R₀ was allowed to vary between 0.001 and 1 mg m⁻² s⁻¹. The same minimum number of observations within a 3 d moving window was used for determining α 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.5 and −0.00001 mg µmol⁻¹ and −5.0 and −0.00001 mg m⁻² s⁻¹, respectively.

Appendix B: Gap-filling of energy fluxes

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 local time, UTC+2). 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 daytime Q_H (R_n>0) were filled with monthly linear regression with R_n. The nighttime gaps in Q_H (R_n<0) were filled with the corresponding R_n values. The gaps in the daytime Q_E were filled in such a way that the monthly mean energy balance closure was achieved. The nighttime gaps in Q_E were set to 0.

Appendix C: Management effect on fluxes

The immediate effect of management on the measured NEE and WUE was investigated by comparing the mean values of 5 d before and after the management day (Table C1).

Table C1 Mean flux and meteorological conditions 5 d before and after management. The management day is not included. NA: not available.