Eddy covariance diverges from inventory and model-based estimates of soil organic carbon stock change across five crop rotations in Belgian cropland

Quantifying long-term soil organic carbon (SOC) dynamics in croplands is critical to evaluate their climate mitigation potential. We assessed SOC stock changes over 21 years at the BE-Lon ICOS cropland (Belgium) by comparing decadal soil inventories, RothC model simulations, and eddy covariance (EC) flux measurements, supported by comprehensive uncertainty propagation. Soil inventories indicated no significant change in SOC stocks, while RothC simulations indicated a gain of 20.4 ± 34.7 g C m⁻² yr⁻¹. In contrast, EC-derived carbon budgets suggested mean annual SOC losses of 69.8 ± 19.2 g C m⁻² yr⁻¹. Uncertainty in the flux-based approach was distributed evenly among spectral corrections, friction velocity thresholds, and lateral carbon flux estimates; yet it remained substantially smaller than the cumulative flux signal. A systematic positive discrepancy between EC-derived budget and RothC simulations suggests an underestimation of net carbon sequestration by the EC technique, which is consistent with the well-documented energy balance non-closure at cropland sites. Despite these methodological discrepancies, both EC and RothC consistently highlighted the critical role of management practices. Notably, both approaches identified the winter wheat–cover crop sequence as a moderate net carbon source in the absence of organic amendments. Conversely, manure application emerged as the primary driver of SOC stock accumulation. More broadly, this study highlights the necessity of combining process-based models, multi-decadal inventories, and EC flux measurements to reliably constrain SOC budgets and identify robust management drivers in cropland ecosystems.

1 Introduction

The management of soil organic carbon (SOC) plays a dual role in sustaining agricultural productivity and regulating the global carbon cycle. While soil represents the largest carbon pool within the terrestrial biosphere, storing two to four times more carbon (C) than the combined total of vegetation and the atmosphere (Lehmann and Kleber, 2015; Terrer et al., 2021), even minor changes in SOC content can either exacerbate or mitigate climate change through their effects on atmospheric CO₂ concentrations (Köchy et al., 2015). As a vital component of agricultural productivity, high SOC stocks improve soil structure, supply nutrients, enhance water retention, and support crop performance (Lal, 2020; Pan et al., 2009). Furthermore, SOC sequestration is increasingly recognized as a viable negative emission method, capable of offsetting one-sixth to one-third of annual net anthropogenic CO₂ emissions (Friedlingstein et al., 2022; Fuss et al., 2018).

Agricultural soils cover an area of more than 1.2 M ha worldwide, spanning across a wide range of climate zone and edaphic conditions (Potapov et al., 2021). The heterogeneity of agricultural landscapes complicates the assessment of SOC changes over space and time (Penko Seidl and Golobič, 2020; Sulman et al., 2020). In particular, SOC dynamics in croplands are influenced by agricultural practices, including tillage, crop rotation, residue management, and fertilization, which collectively influence nutrient cycling, soil fertility, and greenhouse gas emissions (Ren et al., 2020). Comprehensive and robust methods for SOC assessment are essential to support sustainable agricultural management and to mitigate climate change through the implementation of monitoring, reporting, and verification (MRV) tools (Smith et al., 2020).

Traditional SOC assessment methods rely on direct soil sampling and laboratory analyses, providing detailed but spatially and temporally limited estimates of carbon stocks, lacking the temporal resolution required to capture SOC dynamics (Raffeld et al., 2024; Smith et al., 2020). Additionally, the variability of bulk density and SOC concentration on croplands may also complicate the assessment of SOC stock changes from soil inventories (Goidts et al., 2009). Flux tower measurements offer an alternative method for assessing CO₂ exchange at the ecosystem–atmosphere interface through the eddy covariance (EC) technique, which delivers high-frequency observations of net CO₂ fluxes across multiple interacting pools, referred to as net ecosystem exchange (NEE) (Baldocchi, 2014). To date, over 800 site-years of data from more than 130 cropland sites have been collected by the FLUXNET network (Pastorello et al., 2020), thereby offering powerful datasets for studying carbon dynamics in croplands. Notably, net biome productivity (NBP) can be estimated from NEE and net lateral CO₂ fluxes from agricultural activities (e.g., manure application and harvesting) and reflects changes in SOC stocks when dissolved organic and inorganic carbon fluxes, and changes in vegetation C stocks are negligible (Aubinet et al., 2009). However, the extensive post-processing of EC data can introduce additional uncertainties (Aubinet et al., 2012; Hemes et al., 2021). Although substantial progress has been achieved through the development of standardized quality-control procedures in centralized data-processing pipelines (Pastorello et al., 2020), comprehensive uncertainty analysis across all steps of the EC data processing chain remains essential to ensure robust estimates of carbon fluxes and derived budgets. Process-based SOC models simulate soil carbon cycling processes and their interactions with climate, soil properties and management practices to predict SOC dynamics over time under changing environmental and land-use conditions (Smith et al., 2020). The Rothamsted Carbon Model (RothC – Coleman and Jenkinson, 1996) has been widely used under varying land-use, management scenarios, and climate conditions to predict SOC stocks at local (Dechow et al., 2019; Geremew et al., 2024; Ingwersen et al., 2024; Zimmermann et al., 2007) and large scales (Gutierrez et al., 2023; Jordon and Smith, 2022; Morais et al., 2019; Muhammed et al., 2018). RothC is currently used for estimating national-scale CO₂ emissions and removals as a Tier 3 approach in Japan, Switzerland and Australia (Forrester et al., 2024, 2024; Wüst-Galley et al., 2020; Yokozawa et al., 2010), is recommended by the Food and Agriculture Organization of the United Nations (FAO) as a baseline model for mapping SOC sequestration potential (Peralta et al., 2022), and has been integrated into Earth system models (Clark et al., 2011). Yet, predictions of SOC remains strongly sensitive to model structure (Wieder et al., 2014), initial conditions (Todd-Brown et al., 2013) and classical parameterization (Bruni et al., 2022; Hararuk et al., 2014), potentially biasing simulated soil C responses to climate change (Shi et al., 2018). Therefore, uncertainty quantification in soil C models remains critical for assessing the reliability of SOC predictions.

This study assesses changes in SOC stocks over 21 years in a typical 4-year crop rotation at the Lonzée EC station in Wallonia, Belgium. The joint comparison of soil inventories, EC-derived carbon budgets, and RothC simulations provides a unique opportunity to evaluate the consistency of independent approaches, identify potential biases, and improve confidence in SOC stock change estimates. The objectives were to (1) quantify long-term SOC stock changes across methods, (2) assess associated uncertainties, and (3) examine the influence of agricultural management practices on SOC stock dynamics.

2 Material and methods

2.1 Study site

The study site is located in Lonzée, Belgium (50°33^′05.8^′′ N 4°44^′46.4^′′ E, 167 m a.s.l.). The cropland has been cultivated for more than 80 years. Since 2004, sugar beet (SB – Beta vulgaris L.), winter wheat (WW – Triticum aestivum L.) and seed potato (SP – Solanum tuberosum L.) have been grown in a four-year rotation system typical of central Belgium under conventional tillage. Fertilization varied with crop and year but was dominated by repeated nitrogen (N) applications, mainly as liquid N, nitrate, ammonitrate or urea-based fertilizers. The average mineral and organic nitrogen fertilization were respectively 161 kg N_min ha⁻¹ yr⁻¹ and 23.2 kg N_org ha⁻¹ yr⁻¹. In addition, potassium, calcium, sulfur, boron, and magnesium were applied intermittently and some rotations also received organic amendments such as cow manure and lime-based products. Farming activities associated with these rotations are given in Table S1.

The Lonzée site has been part of the FLUXNET network since 2004 and was labelled “class 2” station in the Integrated Carbon Observation System (ICOS) network in 2017 (site code BE-Lon). The climate is temperate oceanic (Köppen classification Cfb) with annual mean precipitation and temperature of 742 mm and 10.2 °C respectively. The soil is classified as a Stagnic Luvisol (WRB classification, Mantel et al., 2023). Texture shows a clear vertical gradient: silt dominates the topsoil (>80 % in 0–30 cm), sand content decreases with depth, and clay content increases correspondingly. SOC content follows the typical exponential decline with depth observed in ploughed agricultural soils, from 10.9 g C kg⁻¹ in 0–30 cm to 3.96 and 1.80 g C kg⁻¹ in the 30–60 and 60–100 cm layers, respectively (measured in 2017; Table A1).

This paper focuses on the 2004–2024 period, thereby including five rotations starting in April 2004. The first (2004–2008), second (2008–2012) and fourth (2016–2020) rotations followed the typical scheme: SB-WW-SP-WW. Cover crops (mustard MU – Sinapis alba L.; oat OAT – Avena sativa L.; faba bean FAB – Vicia faba L., phacelia PHA – Phacelia tanacetifolia Benth., clover CLO – Trifolium repens L.) were sown after harvesting WW except in 2005, 2007 and 2011. The third (2012–2016) and fifth (2020–2024) rotations diverged from this cropping sequence: in 2012 silage maize (MA – Zea mays L.) replaced SB, and spinach (SPI – Spinacia oleracea L.) and green dwarf beans (GWB – Phaseolus vulgaris L. var. nanus) were sown in 2021 to achieve higher economic returns. Due to technical constraints, the spinach crop was not harvested and was instead crushed and incorporated into the soil in May 2021. Consequently, we classified it as a cover crop in our analysis, as all aboveground biomass were returned to the soil. Green dwarf beans were subsequently sown in August 2021 and harvested in October 2021. Consequently, rotation 5 was extended beyond the typical cycle, ending with the crushing of a mix of phacelia and clover as cover crops after winter wheat harvest until the end of 2024. Due to lack of data availability, the three months of bare soil preceding the end of rotation 5 from January 2025 to March 2025 were not included in the analysis. The study period therefore ended with phacelia-clover crushing in December 2024. The crop calendar is given in Fig. 1.

Figure 1 Summary of five crop rotations (2004–2024) at BE-Lon. Each crop period spans from the sowing of one crop to the sowing of the subsequent crop. Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean. Rotations ran from 31 March to 31 March four years later, except Rotation 5, which ended on 31 December 2024. Cover crops included: oat–faba bean (2019–2020), spinach (2021), and phacelia–clover (2024); mustard was used in all other years. FYM is farmyard manure. Hatches show bare soil periods between harvest/crushing and sowing/planting events.

2.2 Data processing

The Lonzée station is equipped with a sonic anemometer (Solent Research R3 before 2014 and Solent Research HS-50 after 2014, Gill Instruments, Lymington, UK) and an infrared gas analyzer (LI-7000 before 2014 and LI-7200 after 2014, LI-COR Inc., Lincoln, NE, US) allowing to quantify the turbulent exchange of CO₂ (NEE), latent heat (LE) and sensible heat (H) from high-frequency (10–20 Hz) measurements of vertical wind speed and gas concentrations. The complete list of sensors is given in Table S2. All variables used in this study, except CO₂ fluxes, can be found in Dumont et al. (2026a). The processing of CO₂ flux data, which is different than the standard FLUXNET pipeline, is described below.

Half-hourly fluxes were calculated from raw data following the methodology described in Aubinet et al. (2012). Block averaging was used for the estimation of fluctuations, double rotation was applied to align the wind vector to the streamlines, and time lag was estimated by covariance maximization. One important step in flux computation consists of applying a spectral correction on CO₂ flux observations. In particular, the EC system acts as a low pass filter due to its limited ability to resolve small-scale variations in vertical wind speed and CO₂ concentration, resulting in attenuation of the measured co-spectra at high frequencies (Aubinet et al., 2012). This signal loss is compensated during post-processing by means of a correction factor (CF), which is used to estimate the true, unattenuated covariance from the measured covariance. Assuming co-spectral similarity among scalar fluxes, CF can be estimated empirically as the ratio of the integrated unbiased wind speed (w_s) – sonic temperature (T_s) co-spectrum to the integrated biased w_s – CO₂ concentration co-spectrum. The latter is modeled as the product of a transfer function and the unbiased w_s–T_s co-spectrum. The transfer function was obtained by fitting a Lorentzian function, characterized by a single fitting parameter corresponding to the cut-off frequency (f_c), to the ratio between high-quality biased and unbiased co-spectra (Ibrom et al., 2007). Beyond this general framework, several implementation choices are possible. In this study, we followed the co-spectral approach described in Faurès et al. (2026). Co-spectra were preferred over spectra for the estimation of f_c. Recent developments by Peltola et al. (2021), highlighting the need to account for an additional filtering linked to the lag search by covariance maximization, were also included. The estimation of f_c and CF was based on binned, class-averaged co-spectra, with atmospheric stability and wind speed used as sorting variables. CF were then used to correct half-hourly CO₂ for compensating high-frequency losses. Co-spectra were computed using the EddyPro software (LI-COR Biosciences, Lincoln, NE, USA, v7.0.9), and subsequent analysis was performed with the Python tool FreqCor v1.3 (see section “Code and data availability”).

Because f_c and CF also depend on the effective measurement height and on instrumental setup characteristics – such as the distance between the sonic anemometer sensing volume and the gas inlet, flow rate, tube length, and filter and rain cap configurations – the analysis was conducted separately for periods during which these characteristics remained constant. The effective measurement height was discretized into three classes split around 1.5 and 2.3 m, resulting in 34 sub-periods over the 21-year study period. Additional information on uncertainty calculations can be found in Sect. 2.4.

After spectral corrections, post-processing of NEE data included (i) determination of friction velocity (u^*) thresholds, and (ii) gap-filling of missing or discarded data under low turbulence, following the REddyProc routine (version 1.3.3; Wutzler et al., 2018). Compared to the typical REddyProc workflow, we implemented several adaptations to account for the crop rotation structure of the Lonzée site. Instead of relying on fixed three-months yearly seasons, we defined a seasonal stratification based on harvest and sowing dates of each crop, reflecting the aerodynamic characteristics of each cropping period. Within each sowing-to-harvest sequence, the following typical data treatment was applied. First, u^* thresholds were determined using the moving point detection method (Papale et al., 2006). Second, to propagate the uncertainty associated with the u^* threshold estimation, five equiprobable u^* scenarios corresponding to the 5th, 16th, 50th, 84th and 95th percentiles of the bootstrapped u^* threshold distribution (n=200 bootstrap samples) were generated. Third, after discarding data under low turbulence (i.e., $u^{*}<u^{*}$ thresholds), gap-filling was performed by applying the standard MDS (Marginal Distribution Sampling) algorithm (Reichstein et al., 2005), using shortwave incoming radiation, air temperature, and vapor pressure deficit as auxiliary meteorological variables. The MDS method consists in replacing missing data with the average of half-hourly CO₂ flux measurements under similar meteorological variables within sliding time-windows. The MDS method was applied for each u^* threshold scenario, yielding an ensemble of five gap-filled NEE time series from which half-hourly uncertainty envelopes (5th–95th percentile range) were derived.

Unlike standard implementations of REddyProc, which treat the full time series as a homogeneous pool for MDS gap-filling, this pipeline segmented the data by crop group prior to gap-filling. This segmentation reduced pool sizes at crop transitions or at the edges of the dataset, necessitating an alternative approach where standard MDS would have failed for gap-filling discarded or missing data. For instance, at the onset of a winter wheat period, only a few days of observations are available within that crop group, making it impossible for MDS to identify a sufficient number of meteorologically similar half-hours within the required window. For these cases, we implemented a fallback that extended the search of similar conditions periods to all years of the dataset while preserving the crop-group consistency. Candidate observations were drawn from the same crop group across the full record and matched using the standard MDS meteorological similarity criteria. The gap-filled value was the arithmetic mean of all retained observations, and a quality flag (fqc = 4–6) was assigned to reflect the level of constraint relaxation, extending the REddyProc flag convention (fqc 0: measured data, fqc 1–3, gap-filled depending on the availability of meteorological drivers). Of the 363 841 half-hourly timesteps retained in the analysis, 67.0 % corresponded to fqc = 0, 29.4 % to fqc = 1–3, and 3.52 % to fqc = 4–6.

2.3 EC-based C budget

Assuming that all dissolved organic and inorganic carbon in precipitation, runoff, percolation and drainage are negligible, the SOC stock budget of an agricultural field over a specific time interval (t to t^′) can be written as (Baker and Griffis, 2005; Zhou et al., 2025):

$$\begin{array}{}\text{(1)}& \mathrm{NBP}=\sum _t^{t^{\prime }}\left(\mathrm{NEE}-{\mathrm{C}}_{\mathrm{imp}}+{\mathrm{C}}_{\exp }\right)\end{array}$$

with NBP the net biome productivity, C_imp the carbon imported (e.g., farmyard manure), C_exp the carbon exported (e.g., harvest), all expressed in mass density (kg C m⁻²). Carbon gain was counted negative and carbon loss positive, following the EC flux sign convention. Since the study period encompasses five full crop rotations, we assumed that residue-derived carbon inputs were effectively balanced across the study period and unlikely to introduce systematic bias in the long-term carbon budget estimates.

The time interval from t to t^′ in Eq. (1) extends from the sowing of one crop to the sowing of the subsequent crop. This sowing-to-sowing definition ensured that most of CO₂ emissions arising from post-harvest crop residue decomposition were attributed to the crop of origin rather than to the following crop, thereby enabling a more robust estimation of the influence of crop species on SOC stock variations. The calculation of the fraction of crop residues that have been decomposed between harvest of the crop and sowing of the next crop will be described in Sect. 2.6.2. Sowing-to-sowing periods also enable a direct comparison between EC- and RothC-based SOC stocks. More information is available in Sect. 2.8.

C_imp from farming practices consisted in the addition of farmyard manure (FYM), limes (SL) and seeds or mother tubers. FYM was added in 2009 (Rotation 2), 2015 (Rotation 3), 2021 and 2024 (Rotation 5, Table 1). C added from FYM (C_imp,FYM) was calculated as follows:

$$\begin{array}{}\text{(2)}& C_{\mathrm{imp},\mathrm{FYM}}={\mathrm{FM}}_{\mathrm{FYM}}\cdot {\left[\mathrm{C}\right]}_{\mathrm{FYM}}\cdot f_{\mathrm{FM}\mathrm{2}\mathrm{DM}}\end{array}$$

where FM_FYM is fresh matter FM of FYM (kg m⁻²), [C]_FYM is C content (kg C kg⁻¹ of dry matter DM) and f_FM2DM is the conversion factor from FM to DM, which is the inverse of the dry matter content (kg DM kg⁻¹). [C]_FYMand f_FM2DM were measured in 2009 and 2021. In practice, samples of FYM were collected at three different depths within the heap of manure (3 repetitions) and dried at 65 °C to determine f_FM2DM. [C]_FYM was measured by dry combustion (ISO 10694:1995). f_FM2DM was 0.29 in 2009 and 0.24 ± 0.07 in 2021 (no standard deviation on f_FM2DM available in 2009). [C]_FYM was 0.33 ± 0.02 and 0.37 ± 0.02 in 2009 and 2021 respectively [C]_FYM and f_FM2DM were therefore set to 0.35 ± 0.01 and to 0.27 ± 0.07 (average of the two campaigns assuming independent measurements). C added from slimes (C_imp,SL) was calculated as follows:

$$\begin{array}{}\text{(3)}& C_{\mathrm{imp},\mathrm{SL}}=m_{\mathrm{SL}}\cdot {\displaystyle \frac{M_{\mathrm{C}}}{M_{\mathrm{SL}}}}\end{array}$$

where m_SL is the mass density of slimes added to the cropland (kg m⁻²), M_C is the molar mass of C, and M_SL is the molar mass of lime (either CaCO₃ or MgCO₃ – kg mol⁻¹). C imported from sowing (C_imp,SE) was calculated as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{C}}_{\mathrm{imp},\mathrm{SE}}=m_{\mathrm{SE}}\cdot \mathrm{SD}\cdot {\left[\mathrm{C}\right]}_{\mathrm{SE}}\end{array}$$

with SD the density of seeds (m⁻²), m_SE the mass of seeds (kg) and [C]_SE the carbon content of seeds (kg C kg⁻¹ DM). m_SE was set to 20 mg (Heimbach, 2018) and [C]_SE was set to 45 % (Ma et al., 2018). Exports of C (C_exp) were calculated as follows:

$$\begin{array}{}\text{(5)}& {\mathrm{C}}_{\exp }={\mathrm{DM}}_{\mathrm{H}}\cdot {\left[\mathrm{C}\right]}_{\mathrm{H}}\end{array}$$

with DM_H the dry matter of harvested material (kg m⁻²) and [C]_H the C content (kg C kg⁻¹ DM) of harvested material (being tubers for SP and SB, grains and a variable fraction of straws for WW, grains and straws for MA, grain legumes for GDB). DM_H was obtained after drying fresh matter harvested in an oven at 65 °C until constant weight. [C]_H was measured by dry combustion (ISO 10694:1995). In 2007, 2013 and 2017, missing harvested straw used to calculate DM_H of winter wheat was estimated from a linear regression based on total dry biomass, fitted across all available observations for winter wheat (R² = 0.91, data not shown). DM_H was then calculated as the sum of predicted harvested straw and measured harvested grain dry matter. C_imp,SEof potato was calculated from dry matter of mother tubers and C content similarly to Eq. (5). Missing measurements of [C]_H were replaced by the average of measurements for similar crops. C_imp associated to FYM and SL are shown in Table 1. DM_H, [C]_H and SD used to calculate C_imp and C_exp are provided in Table A2. FM_FYM, m_SL and SD were communicated by the farmer. All other ancillary data collected at BE-Lon follows ICOS standard protocols (Gielen et al., 2018) and can be found in Dumont et al. (2026b) for the period 2017–2024. For the period 2004–2017, the protocol for ancillary data collection is described in Buysse et al. (2017).

Table 1 Carbon inputs added over the five-rotation period at BE-Lon. FM_FYM is fresh matter of farmyard manure (FYM), m_SL is the density of slimes (SL) and C_imp is the corresponding input of carbon (Eqs. 2–4). All uncertainties are cumulative standard deviation.

2.4 Uncertainty on EC-based C budget

The half-hourly uncertainty on NBP (σ_NBP,30) was calculated as the square root of the sum of squares of all individual half-hourly uncertainties on NEE (σ_NEE,30), C_imp (${\mathit{\sigma }}_{{\mathrm{C}}_{\mathrm{imp}},\mathrm{30}}$) and C_exp (${\mathit{\sigma }}_{{\mathrm{C}}_{\exp },\mathrm{30}}$) as:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{NBP},\mathrm{30}}=\sqrt{{\left({\mathit{\sigma }}_{\mathrm{NEE},\mathrm{30}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{{\mathrm{C}}_{\mathrm{imp}},\mathrm{30}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{{\mathrm{C}}_{\exp },\mathrm{30}}\right)}^{\mathrm{2}}}\end{array}$$

Equation (6) explicitly considers that the covariance between all uncertainties is negligible, that is, errors are assumed independent. σ_NEE,30 was calculated as the combination of uncertainties related to the low turbulence filtering procedure (${\mathit{\sigma }}_{\mathrm{NEE},u^{*},\mathrm{30}}$), to the spectral correction procedure (${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},\mathrm{30}}$) and to the random error in NEE measurement (${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{meas},\mathrm{30}}$) (Aubinet et al., 2012). As the MDS approach generally yields results comparable to alternative methods (Teng et al., 2020; Zhao and Huang, 2015), particularly at sites below 60° N and in the absence of prolonged data gaps (Vekuri et al., 2023, 2025), the uncertainty in gap-filling was neglected in the analysis. Equation (6) then becomes:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{NBP},\mathrm{30}}=\sqrt{\begin{array}{c}{\left({\mathit{\sigma }}_{\mathrm{NEE},\mathrm{meas},\mathrm{30}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{\mathrm{NEE},u^{*},\mathrm{30}}\right)}^{\mathrm{2}}\\ +{\left({\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},\mathrm{30}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{{\mathrm{C}}_{\mathrm{imp}},\mathrm{30}}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{{\mathrm{C}}_{\exp },\mathrm{30}}\right)}^{\mathrm{2}}\end{array}}\end{array}$$

The calculation of each term of Eq. (7) and the propagation of uncertainty over long-term budgets is described below.

2.4.1 Uncertainty on lateral C fluxes

${\mathit{\sigma }}_{{\mathrm{C}}_{\mathrm{imp}},\mathrm{30}}$ and ${\mathit{\sigma }}_{{\mathrm{C}}_{\exp },\mathrm{30}}$ were calculated using Eqs. (2)–(5) and the Gaussian error propagation (Ku, 1965), giving for C_exp as an example:

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{{\mathrm{C}}_{\exp },\mathrm{30}}=\sqrt{{\left({\left[\mathrm{C}\right]}_{\mathrm{H}}\cdot {\mathit{\sigma }}_{{\mathrm{DM}}_{\mathrm{H}}}\right)}^{\mathrm{2}}+{\left({\mathrm{DM}}_{\mathrm{H}}\cdot {\mathit{\sigma }}_{{\left[\mathrm{C}\right]}_{\mathrm{H}}}\right)}^{\mathrm{2}}}\end{array}$$

For ${\mathit{\sigma }}_{{\mathrm{C}}_{\mathrm{imp}},\mathrm{30}},$no standard deviations on FM_FYM (Eq. 2) and m_SL (Eq. 3) were available. The uncertainty on FM_FYM was set to ±10 % to account for optimistic typical field-scale deviations caused by variations in manure density, volume estimation, and spreader calibration (Thirion and Chabot, 2003). These uncertainties were attributed to the median half-hour of the day (i.e., 12 h GMT+1) when measurements were performed.

2.4.2 Uncertainty on NEE

${\mathit{\sigma }}_{\mathrm{NEE},u^{*},\mathrm{30}}$ was calculated as follows (Pastorello et al., 2020):

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_{\mathrm{NEE},u^{*},\mathrm{30}}={\displaystyle \frac{\mathrm{NEE}\mathrm{\_}\mathrm{U}\mathrm{84}\mathrm{\_}\mathrm{f}-\mathrm{NEE}\mathrm{\_}\mathrm{U}\mathrm{16}\mathrm{\_}\mathrm{f}}{\mathrm{2}}}\end{array}$$

where NEE_U84_f and NEE_U16_f are respectively the half-hourly gap-filled NEE after applying a u^* threshold corresponding to the 84th and 16th percentile of the threshold distribution within a sowing-to-harvest period (see Sect. 2.2). ${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{meas},\mathrm{30}}$ was calculated as the standard deviation of the measured NEE fluxes following the method described in Hollinger and Richardson (2005).

${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},\mathrm{30}}$ corresponds to the total uncertainty related to the spectral correction applied on CO₂ flux measurements, which consists of two components: (i) one related to the cut-off frequency resulting from the fitting of the Lorentzian function, referred to as the transfer function uncertainty σ_tf, and (ii) one related to the aggregation of correction factors within wind speed and stability classes, referred to as the spectral class uncertainty σ_cl. σ_tf corresponds to the square root of the diagonal of the covariance matrix associated with the fitted Lorentzian parameters, while σ_cl was calculated as the standard deviation of the population of correction factors CF for a given wind speed and stability class within a spectral correction period. These two uncertainty sources are physically independent: σ_tf depends solely on the uncertainty in the cut-off frequency f_c derived from the Lorentzian fitting procedure, while σ_cl reflects the residual variability of individual correction factors within each class, which can be attributed to other factors such as like crop height or wind direction.

2.4.3 Propagation of uncertainties

The general framework of uncertainty propagation posits that systematic errors accumulate linearly and random errors accumulate in quadrature (Richardson et al., 2012). Therefore, the cumulative uncertainty on NEE and ultimately on NBP strongly depends on the categorization of each uncertainty term as systematic or random.

For ${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},\mathrm{30}}$, we considered each spectral correction period as an independent correction unit, within which we applied a nested error structure reflecting the nature of the spectral correction procedure. Within each wind speed – stability class combination (k) for a given spectral correction period (p), since σ_tf,k originates from a single cut-off frequency estimate derived from the Lorentzian fitting procedure and applied identically to all half-hours of class k, it was treated as systematic and propagated linearly over the cumulated NEE of the class (Richardson et al., 2012). The uncertainty on NEE related to the cut-off frequency for each class (${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{tf},k}$) was therefore calculated as:

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{NEE},\mathrm{tf},k}={\mathit{\sigma }}_{\mathrm{tf},k}\cdot \sum _{t\in k}\left|{\mathrm{NEE}}_t\right|\end{array}$$

In contrast, ${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{cl},k}$ was treated as random, assuming that variability of co-spectra within a wind speed – stability class is uncorrelated between half-hours, and was therefore propagated quadratically:

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_{\mathrm{NEE},\mathrm{cl},k}={\mathit{\sigma }}_{\mathrm{cl},k}\cdot \sqrt{\sum _{t\in k}{\left({\mathrm{NEE}}_t\right)}^{\mathrm{2}}}\end{array}$$

The composite uncertainty for each class k was then obtained by combining these two independent contributions quadratically:

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},k}=\sqrt{{\left({\mathit{\sigma }}_{\mathrm{NEE},\mathrm{cl},k}\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{\mathrm{NEE},\mathrm{tf},k}\right)}^{\mathrm{2}}}\end{array}$$

Since wind speed – stability classes capture distinct turbulent regimes, their associated errors can be considered independent across classes. The total standard deviation on cumulated NEE within a spectral correction period p containing n_k wind speed classes was therefore obtained as the square root of the sum of squares of class-level contributions:

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},p}=\sqrt{\sum _{k=\mathrm{1}}^{n_k}{\left({\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},k}\right)}^{\mathrm{2}}}\end{array}$$

For simplicity, the class-level uncertainty σ_NEE,k was redistributed to each half-hour in proportion to its absolute NEE value, normalized by the root-sum-of-squares of NEE within the class. This provides an estimate of the half-hourly uncertainty in NEE due to the spectral correction procedure (${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},\mathrm{30}}$). This redistribution enables that the quadratic sum of half-hourly uncertainties exactly recovers ${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},k}$, and by extension ${\mathit{\sigma }}_{\mathrm{NEE},\mathrm{spcorr},p}$ at every aggregation level.

A very similar approach was applied to the uncertainty on u^* threshold determination. Since a single u^* threshold is applied uniformly to all half-hours within a given crop season, this uncertainty is systematic in nature within each season. Accordingly, errors were accumulated linearly within each season to yield a seasonal u^* uncertainty. Since crop seasons are determined independently and their respective thresholds are estimated from non-overlapping datasets, seasonal uncertainties are treated as independent and combined in quadrature across seasons. Half-hourly values were then redistributed from the seasonal budget in proportion to their individual contribution, preserving quadratic consistency with the total uncertainty related to u^* threshold procedure, and giving ${\mathit{\sigma }}_{\mathrm{NEE},u^{*},\mathrm{30}}$. σ_NBP,30 calculated using Eq. (7) was then propagated in quadrature over a specified period (d) following:

$$\begin{array}{}\text{(14)}& \begin{array}{l}{\mathit{\sigma }}_{\mathrm{NBP},d}=\sqrt{\sum _{t\in d}{\left({\mathit{\sigma }}_{\mathrm{NBP},\mathrm{30},t}\right)}^{\mathrm{2}}}\end{array}\end{array}$$

where t is the number of half-hourly data within the corresponding time window d. As the covariance was neglected, the reported NBP uncertainty should therefore be considered as a higher bound on the true uncertainty.

2.5 Soil inventories

Soil samples were collected in 2007 and 2017. In 2007, the sampling scheme corresponding to a systematic grid covered a rectangular area of 0.6 ha around the EC flux tower (Rodeghiero et al., 2010; Schrumpf et al., 2011). The 61 cores extracted at the node of the sampling grid were each divided into 2–3 pedogenetic horizons, reaching a mean maximum depth of 60 cm. In 2017, the sampling scheme covered the entire cropland area (14 ha), which corresponded to a stratified random sampling where the sampling domain was divided into six strata of equal area. 20 sample locations were selected from the sampling domain from which 100 samples were extracted over 5 depths (0–5, 5–15, 15–30, 30–60, 60–100 cm) following ICOS protocols (Arrouays et al., 2018). As the carbonate content was below 0.70 g C inorganic kg⁻¹ (0.34, 0.36, 0.47 g C inorganic kg⁻¹ over 0–5, 5–15 and 15–30 cm, Table A1), SOC concentration was calculated as the difference between total and carbonate concentration (Arrouays et al., 2018). In 2007, no information on carbonate concentration was available. Since the soil is classified as Stagnic Luvisoil (Mantel et al., 2023) developed under leaching conditions and inherently non-calcareous; carbonate contents are therefore negligible, and total carbon is taken as equivalent to SOC for the 2007 campaign.

SOC stocks were computed on an equivalent soil mass basis (ESM; Lee et al., 2009) using the simpleESM R script (Ferchaud et al., 2023) with reference soil masses of 384 and 436 Mg ha⁻¹ in ∼ 0–30 and ∼ 30–60 cm layers. The mean and standard deviation of SOC stocks over these layers were then computed using a design-based approach (Loustau et al., 2017). Because the soil was relatively homogeneous, a reduced dataset was reconstructed by selecting spatially co-located sampling points from 2007 and 2017 to ensure a reliable comparison between the two datasets (Fig. A1) (Loubet et al., 2026). Differences in SOC stocks over the ∼ 0–30 and ∼ 30–60 cm layers between 2007 and 2017 within the reduced area were quantified using Hedges' g effect size estimates, allowing to correct standardized mean difference for small-sample bias (Hedges, 1981). Confidence intervals (CI) on Hedges' g were obtained using 10 000 nonparametric bootstrap resamples and estimated with the bias-corrected and accelerated method (Kirby and Gerlanc, 2013). More information on the statistical analysis on SOC stock comparison can be found in Loubet et al. (2026).

2.6 Process-based modelling of SOC

2.6.1 The RothC model

The RothC model (Coleman and Jenkinson, 1996) was used to simulate SOC stock changes in the topsoil (i.e., over 0–23 cm, as the model was calibrated and validated in this soil layer) between 2004 and 2024. RothC divides SOC into five conceptual pools: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO), humified organic matter (HUM) and inert organic matter (IOM) pools. SOC stock depletion within each pool are influenced by rate modifying factors which depend on climate conditions, land use and pool characteristics (Coleman and Jenkinson, 1996). RothC requires the following inputs: initial SOC stock, precipitation, open-pan evaporation, air temperature, clay content, decomposability of incoming plant material (expressed as the $\mathrm{DPM}/\mathrm{RPM}$ ratio and set to 1.44 for croplands), soil cover (either 0.6 for vegetated soil or 1.0 for bare soil) and C inputs from farmyard manure and plant biomass (in kg C m⁻²). The latter includes all plant material returned during the cultivation period, including both above ground biomass as crop residues and belowground biomass as C exudates or unharvested root biomass. In practice, fractions of carbon returned to the soil as belowground biomass (C_BGB) were estimated from the sum of plant carbon in root tissues (C_R) and root exudates (C_E) (Bolinder et al., 2007):

$$\begin{array}{}\text{(15)}& {\mathrm{C}}_{\mathrm{BGB}}={\mathrm{C}}_{\mathrm{R}}+{\mathrm{C}}_{\mathrm{E}}\end{array}$$

where C_E represents 65 % of C_R (Bolinder et al., 2007). For silage maize and winter wheat, C_R was estimated from carbon exported from the main agricultural product (C_P), shoot-to-root ratio (S:R) and harvest index (HI) (Bolinder et al., 2007; Keel et al., 2017):

$$\begin{array}{}\text{(16)}& {\mathrm{C}}_{\mathrm{R}}={\displaystyle \frac{C_{\mathrm{P}}}{S:R\cdot \mathrm{HI}}}\end{array}$$

HI was calculated from the ratio of harvested grain to total biomass harvested and S:R was quantified during a 2022 pre-harvest field campaign on winter wheat. Measurements were conducted at four randomly selected locations across the cropland, following established protocols for both shoot biomass (Gielen et al., 2018) and root biomass (Frasier et al., 2016; Noordwijk et al., 1985). Root biomass was measured by collecting soil samples at six depth intervals within the 0–100 cm layer using hand augers. Roots were separated from soil using a saline solution (NaCl), followed by sieving (2, 1 and 0.5 mm). The samples were then dried at 65 °C until constant weight was achieved. Total root biomass was calculated by normalizing root weight by the soil sampling area across the entire soil profile. Shoot biomass was determined by harvesting 5 randomly selected plants at each sampling location. All aboveground organs were dried at 65 °C until reaching constant weight. Final shoot biomass was calculated as the ratio of dry organ weight to the sampled soil surface area. S:R was 12.3 which in the upper end of the values reported by Bolinder et al. (2007). S:R for silage maize was assumed to be equivalent to that of winter wheat  (Bolinder et al., 2007). For the other crops, Bolinder coefficients were used to determine C_R and C_E from yield or from crushed aboveground biomass for cover crops (Bolinder et al., 2007). C_BGB was rescaled to 23 cm to match RothC layers definition (Dechow et al., 2019; Keel et al., 2017).

2.6.2 Initialization and simulation phases

Reactive carbon pools in RothC were initialized as a function of clay content and initial SOC stock (Falloon et al., 1998; Morais et al., 2019; Weihermüller et al., 2013). The initial SOC stock in 2004 was set to the 2007 inventory measurement, since significant changes in SOC stocks are generally detectable only after at least 10 years (Schrumpf et al., 2011; Smith, 2004). All model parameters were chosen accordingly to the original version of RothC (Coleman and Jenkinson, 1996). The decomposition rate constants related to pool characteristics were not modified (10 yr⁻¹ for DPM, 0.3 yr⁻¹ for RPM, 0.66 yr⁻¹ for BIO and 0.02 yr⁻¹ for HUM). The model was then run between 2004 and 2024 using observed meteorological data, C inputs from plant biomass and from FYM. Potential evapotranspiration was calculated using the FAO Penman-Monteith equation (Pereira et al., 2025; Cai et al., 2007) based on wind speed, vapor pressure deficit, and net radiation. Open-pan evapotranspiration was calculated by dividing potential evapotranspiration by 0.75 (Coleman and Jenkinson, 1996). For winter wheat, maize and cover crops, carbon inputs were calculated as the sum of measured aboveground crop residues after harvest (Table A3) and estimated belowground biomass (C_BGB – Eq. 15). For tuber crops, measured harvest residues were selected as carbon inputs. For each simulated month, the difference between the initial SOC stock was calculated (ΔSOC_RothC).

To draw meaningful conclusions regarding the influence of crop species on SOC stock variations, we quantified the potential carry-over C residues over the sowing-to-sowing crop period by calculating the fraction of the C input allocated to the DPM and RPM pools in RothC (f_res) as follows:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{f}_{\mathrm{res}}& =\mathrm{1}-{\displaystyle \frac{{\mathrm{C}}_{\mathrm{DPM}}\left(t^{\prime }\right)+{\mathrm{C}}_{\mathrm{RPM}}\left(t^{\prime }\right)}{{\mathrm{C}}_{\mathrm{DPM}}\left(t\right)+{\mathrm{C}}_{\mathrm{RPM}}\left(t\right)}}\\ & =\mathrm{1}-\left(p_{\mathrm{DPM}}\cdot e^{-k_{\mathrm{DPM}}{\sum }_{i=t}^{t^{\prime }}{\mathit{\xi }}_i}+p_{\mathrm{RPM}}\cdot e^{-k_{\mathrm{RPM}}{\sum }_{i=t}^{t^{\prime }}{\mathit{\xi }}_i}\right)\end{array}\end{array}$$

where C_DPM and C_RPM are the C in the DPM and RPM pools respectively, ξ_i is the product of RothC rate modifiers for temperature, soil moisture and soil cover over a month of index i between t and t^′, k_DPM and k_RPMare the monthly DPM and RPM decomposition rate constants (0.83 and 0.025 per month respectively), and $p_{\mathrm{DPM}}=\frac{\mathrm{DR}}{\mathrm{DR}+\mathrm{1}}$ and $p_{\mathrm{RPM}}=\frac{\mathrm{1}}{\mathrm{DR}+\mathrm{1}}$ are the fraction of C residue attributed to the DPM and RPM pool respectively with DR the DPM:RPM ratio of 1.44. The summation is taken over the months included in the sowing-to-sowing period from t to t^′. We focused on DPM and RPM as it corresponds to the pools most likely to generate a short-term carry-over of newly added residue C across successive sowing-to-sowing periods. Since BIO and HUM have substantially lower decomposition rate constants and represent slower C dynamics, they were not used in the calculation of f_res. Periods with low f_res were interpreted as periods with substantial carry-over of DPM and RPM C beyond the end of the crop period. In these cases, the attribution of residue-related C turnover and associated C losses to an individual crop period may be biased, and the resulting sink or source classification should be interpreted with caution.

2.7 Uncertainty on RothC simulations

Monte Carlo (MC) simulations represents a powerful probabilistic framework for quantifying the uncertainties related to inputs, parameters, model structure, and initial conditions (Pierson et al., 2022). MC methodology was applied to RothC by propagating uncertainties through the model via repeated random sampling in probability distribution over 10 000 simulations (Morais et al., 2019; Stamati et al., 2013). The process notably involves selecting the range and probability distributions of selected factors, generating random samples, running the model for each sample, and then conducting uncertainty analysis (Saltelli et al., 2000). Model-input uncertainties were propagated by independently sampling monthly values from Gaussian distributions centred on their measured values. For air temperature and precipitation, standard deviations were derived from sensor accuracy (Vaisala, 2025; MicroStep-MIS, 2026). Monthly plant-derived C inputs, FYM C inputs, and initial SOC stock were sampled using their measurement-based standard deviations. Uncertainty in open-pan evapotranspiration and clay content was neglected. Finally, statistical distributions of RothC parameters (decomposition rates, temperature and bare soil moisture deficit rate modifiers) were also defined based on a Gaussian distribution but within a 20 % margin, which is a typical approach when uncertainty around parameters is not known (Delhez et al., 2025). For each MC simulation, the monthly total SOC and the corresponding difference with initial SOC were calculated. For each month, mean, median, and standard deviation, were calculated across the ensemble of Monte Carlo simulations, yielding a distinct uncertainty estimate for every monthly time step. RothC was run in R (R Core Team, 2023) using the SoilR package (version 1.2; Sierra et al., 2012).

2.8 Comparison between cumulative NBP and SOC stocks simulated by RothC

Unlike inventories, which measure changes in SOC stock between two discrete sampling campaigns, cumulative NBP and ΔSOC_RothC provide continuous estimates of SOC stock variations. However, the two methods rely on fundamentally different principles: NBP integrates NEE, capturing the net CO₂ exchange between vegetation and the atmosphere, and thus directly reflecting the net ecosystem carbon balance after accounting for autotrophic and heterotrophic respiration, as well as net lateral C fluxes. In contrast, RothC does not simulate plant uptake explicitly and contains no living-biomass compartment, no plant growth sub-model, and resolves only the decomposition and turnover of organic matter within the soil from prescribed carbon inputs.

This structural difference defines the conditions under which the two quantities can be compared. Over any given interval, the ecosystem carbon balance can be decomposed as the sum of the net changes in soil C stocks and in living biomass C stocks. Since RothC predicts changes in SOC stocks only, EC-based C budgets and RothC simulations can only be compared when changes in living biomass C stocks are negligible over the selected period. We therefore restricted the comparison to sowing-to-sowing periods. At the beginning of the period, no vegetation is present on the field and living biomass C stocks are zero. At the end of the period (sowing of the subsequent crop, which can occur a few months after harvest), living biomass C stocks have also returned to near zero, since crop residues left in the field after harvest enter the soil through ploughing. This step is explicitly considered in RothC, where aboveground residues left on the field are prescribed to the DPM and RPM pools at the month of harvest. Belowground crop residues, consisting of root material built up during the cultivation period, are captured by NBP through NEE and are likewise prescribed to RothC through allometric equations based on yield. Any respiratory flux resulting from harvest residue mineralization is accounted for by both approaches. Carbon losses through DOC and DIC leaching or erosion were assumed negligible at this site and are not accounted for in either method. Overall, the comparison between cumulative NBP and ΔSOC_RothC consisted in two parts: first comparing their dynamics to identify systematic trends or biases over the five-rotation period, and second, the calculation of cumulative NBP and ΔSOC_RothC over individual sowing-to-sowing periods. In particular, the net change in each variable was computed as the difference between cumulative values at the end and beginning of each sowing-to-sowing period.

As RothC operates at a monthly time step while EC measurements are provided at half-hourly resolution, NBP and its associated uncertainty were aggregated to a monthly basis prior to comparison. For each sowing-to-sowing period, cumulative NBP was extracted at the exact day of transition between consecutive periods within the month, thereby capturing the net carbon exchange integrated over the sowing-to-sowing periods, including post-harvest residue mineralization. Because both cumulative NBP and ΔSOC_RothC time series carry correlated uncertainties (each monthly value inheriting the accumulated error of all preceding timesteps) the variance of the difference between the two stocks $\mathrm{\Delta }S=S\left(t\right)-S\left(t_{\mathrm{0}}\right)$ reduces to ${\mathit{\sigma }}^{\mathrm{2}}\left(\mathrm{\Delta }S\right)={\mathit{\sigma }}^{\mathrm{2}}\left(S_t\right)-{\mathit{\sigma }}^{\mathrm{2}}\left(S_{t_{\mathrm{0}}}\right)$. The resulting standard deviation σ(ΔS) was used as the uncertainty estimate for each period. The sum of the difference between cumulative NBP and ΔSOC_RothC variations was calculated within each period to quantify agreement between the two estimates. This comparison helps determine whether both methods converge on the same source–sink pattern for a given crop period. RothC simulations and soil inventories followed the same sign convention as NBP, in which carbon gains were considered negative and carbon losses positive.

3 Results

3.1 Eddy covariance-based carbon budget

Winter wheat was the crop with the longest cumulative cultivation period over the five-rotation period for a total of 129.1 months, followed by cover crop (53.8 months), sugar beet (32.2 months), potato (28.8 months), maize (5.46 months) and green dwarf bean (3.13 months). Carbon exports were highest for winter wheat (5.49 ± 0.13 kg C m⁻²), followed by sugar beet (3.64 ± 0.22 kg C m⁻²), seed potato (1.55 ± 0.01 kg C m⁻²), maize (0.78 ± 0.01 kg C m⁻²), and green dwarf bean (0.04 ± 0.00 kg C m⁻²) (Table 2, Fig. 2).

Table 2 Cumulative values of net ecosystem exchange (NEE), imports of carbon (C_imp), exports of carbon at harvest (C_exp), net biome productivity (NBP) for each rotation and each crop over the five-rotation period. Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean. Crops are mentioned with the corresponding year. All uncertainties are cumulative standard deviation.

Figure 2 Cumulative net ecosystem exchange (NEE), exports of carbon (C_exp), imports of carbon (C_imp) and net biome productivity (NBP) for each crop over the five-rotation period. The error bars on the bar plots show cumulative standard deviations on each budget term. Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean.

Farmyard manure was the main external carbon input across the study period, with applications restricted to Rotations 2, 3 and 5 (Table 2). Cumulative NEE differed markedly among crops: it was strongly negative (i.e. net C uptake) under winter wheat and sugar beet, near zero under potato, silage maize and green dwarf bean, and positive (i.e. net C release) under the cover crop (Table 2). As manure application occurred after winter wheat harvest in 2009, 2015 and 2024 and after cover crop crushing in 2021 (Fig. 1, Table 1), the cumulative NBP of these crops was therefore reduced (Table A3). In particular, the cumulative NBP of winter wheat was negative (−1.02 ± 0.25 kg C m⁻²) while other crops exhibited cumulative NBP values ranging from near equilibrium to a moderate carbon source (Fig. 2, Table 2). Positive cumulative NBP values were observed for maize (0.16 ± 0.05 kg C m⁻²), potato (0.91 ± 0.09 kg C m⁻²), sugar beet (0.32 ± 0.24 kg C m⁻²) and cover crops (1.10 ± 0.16 kg C m⁻²) whereas cumulative NBP of green dwarf bean (−0.02 ± 0.04 kg C m⁻²) did not differ significantly from zero (Table 2).

Over the five rotations (252 months), 11.5 ± 0.25 kg C m⁻² were exported at harvest while −7.72 ± 0.22 kg C m⁻² was sequestrated by vegetation and −2.28 ± 0.21 kg C m⁻² was added as FYM, slimes, seeds and mother tubers, giving a total cumulative NBP of 1.45 ± 0.40 kg C m⁻² (Tables 1, 2). It corresponds to a yearly loss of 69.8 ± 19.2 g C m⁻² on average. The uncertainty in cumulative NBP was distributed relatively evenly across all components of the carbon budget. The cumulative ${\mathit{\sigma }}_{{\mathrm{C}}_{\mathrm{imp}}}$ and ${\mathit{\sigma }}_{{\mathrm{C}}_{\exp }}$ were 0.25 and 0.21 kg C m⁻², respectively. ${\mathit{\sigma }}_{\mathrm{NEE},u^{*}}$ reached 0.17 kg C m⁻² while σ_NEE,spcorr amounted to 0.13 kg C m⁻², leading to σ_NEE of 0.22 kg C m⁻² (Table 2). NBP was the highest for Rotation 4, which is explained by the highest exports of C in the five-rotation period (Fig. 3, Table 2). Despite no cover crops being grown and no manure being applied, Rotation 1 had a cumulative NBP which remained close to equilibrium (Fig. 3, Table 2). Compared to Rotation 1, Rotation 2 received an additional 0.30 kg C m⁻² of C inputs via manure (Fig. 3, Table 2). As explained in Sect. 2.1, three months of bare soil following cover crop crushing were not included in Rotation 5 between January and April 2025, which would have slightly increased the cumulative NBP due to heterotrophic soil respiration. Other rotations showed higher positive cumulative NBP trends despite manure applications in 2009 and 2015 (Fig. 3). Overall, the persistent positive trend across the observation period indicates a long-term carbon source pattern. Cumulative carbon fluxes for each crop-period are given in Table A3.

Figure 3 Cumulative net ecosystem exchange (NEE), exports of carbon (C_exp), imports of carbon (C_imp) and net biome productivity (NBP) for each rotation over the five-rotation period. The shaded curve represents the cumulative standard deviation on NBP (σ_NBP) over each rotation. The error bars on the bar plots show cumulative standard deviations on NEE, C_imp and C_exp. Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean.

3.2 Comparison between RothC simulations, cumulative NBP and inventories

The fraction of mineralized DPM and RPM pools (f_res) calculated from RothC SOC stock simulations, varied widely across crop-periods (Table 3). As C inputs in RothC are attributed to the harvest month and the three preceding months, f_res was driven primarily by the length of the harvest-to-sowing period and the corresponding meteorological conditions (Table 3).

Table 3 Harvest-to-sowing intervals and fraction of mineralization of the DPM and RPM pools (f_res). Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean.

Cover crop periods consistently showed f_res of ∼ 40 %–50 %, as residues were crushed in winter under conditions favourable to early-spring mineralization, except for 2024 due to the truncation of the study period on 31 December. Similarly, winter wheat periods without cover crops (2004–2006, 2006–2008, 2010–2012) yielded harvest-to-sowing interval exceeding 8 months and f_res∼60 %. Conversely, short harvest-to-sowing interval associated with winter wheat (2008–2009, 2014–2015, 2018–2019, 2023–2024), maize (2012), sugar beet (2016), seed potato (2018) and green dwarf bean (2021) resulted in markedly lower mineralized fractions (Table 3). Interestingly, when the short winter wheat harvest-to-sowing intervals of 2008–2009, 2014–2015, 2018–2019 and 2023–2024 were combined with their immediately following cover crop periods to form full rotation sequences, the f_res increased substantially, reaching 62.9 %, 64.6 %, 62.6 %, and 66.6 %, respectively. Taken together, these results demonstrate that the WW–CC sequence consistently achieves a high f_res, strengthening the reliability of the NBP and RothC simulated carbon balance, as the sink or source status of the crop-period is not heavily confounded by undecomposed DPM and RPM pools across interval boundaries. Overall, f_res was strongly and positively correlated with harvest-to-sowing interval length across all crop-periods (Pearson r=0.74, p<0.001), confirming that interval duration is the primary control on DPM mineralization at BE-Lon. No significant additional effect of harvest timing was detected after controlling for interval length (partial $r=-\mathrm{0.16}$, p=0.42).

RothC simulations yielded a final ΔSOC of −0.43 ± 0.73 kg C m⁻² over the five-rotation period (Fig. 4). The large standard deviation reflects an uncertain carbon balance, ranging from moderate sink to source. Increase in ΔSOC_RothC stocks occurred either due to manure application or during the three months preceding and including the harvest month due to crop residues and root exudates (Fig. 4). On the other hand, positive ΔSOC_RothC (i.e., decrease in SOC stocks) notably occurred after harvest, where CO₂ emissions were solely driven by heterotrophic respiration resulting from crop residues (Fig. 4).

Figure 4 Cumulative net biome productivity (NBP) and SOC stock variation simulated by RothC (ΔSOC_RothC). Vertical lines indicate farmyard manure (FYM) application. The shaded curves represent the cumulative standard deviation over the study period. The black circle indicates the difference between the SOC stocks measured in 2007 and 2017 (ΔSOC_inv). Note that the grey curve corresponds to the cumulated NBP per month, which leads to a monthly aggregation of events (e.g., FYM and harvest if occurring the same month).

Cumulative SOC stocks over the ∼ 0–60 cm profile (820 kg m⁻² of fine earth), derived from inventories on an equivalent soil mass basis, were 5.99 ± 0.14 kg C m⁻² in 2007 and 5.93 ± 0.22 kg C m⁻² in 2017. Corresponding stocks in the topsoil (∼ 0–30 cm) were 3.81 ± 0.38 and 3.79 ± 0.21 kg C m⁻², respectively, while subsoil (∼ 30–60 cm) stocks were 2.18 ± 0.48 and 2.17 ± 0.28 kg C m⁻², respectively. The magnitude of change between the two inventories (ΔSOC_inv) was negligible, as indicated by Hedges' g values overlapping zero for both layers (topsoil: −0.05, with a 95 % CI from −0.45 to 0.40; subsoil: −0.02, with a 95 % CI from −0.38 to 0.43). Therefore, ΔSOC_inv was not significantly different from 0 over the period 2007–2017. ΔSOC_RothC in 2017 (−0.28 ± 0.48 kg C m⁻²) was of similar magnitude to the inventory change, whereas cumulative NBP over the same period was higher than 1.0 kg C m⁻² (Fig. 4).

By the end of the five-rotation period, the discrepancy between ΔSOC_RothC and cumulative NBP reached 1.88 ± 0.83 kg C m⁻². The comparison between variations of cumulative NBP and ΔSOC_RothC for each crop periods showed a systematic positive discrepancy between the two variables as highlighted by the cumulative difference for maize, seed potato and winter wheat with respectively 0.25 ± 0.08, 0.63 ± 0.08, and 1.55 ± 0.54 kg C m⁻² (Fig. 5). For sugar beet, the cumulative difference was 0.28 ± 0.35 kg C m⁻², which was strongly influenced by the year 2008 with a negative cumulative NBP of −0.28 ± 0.35 kg C m⁻² while ΔSOC_RothC was close to zero (Fig. 5, Table A3). In particular, the cumulative NEE for this crop over the active vegetation period in 2008 (from sowing to harvest) was −1.26 ± 0.03 kg C m⁻², which was 36 %, 85 %, and 129 % higher in absolute value than in 2004 (−0.81 ± 0.05 kg C m⁻²), 2016 (−0.68 ± 0.02 kg C m⁻²), and 2021 (−0.55 ± 0.08 kg C m⁻²), respectively (data not shown). In the meantime, the yield and crop residues did not follow the same magnitude (Table A3), which might indicate an underestimation of biomass measured at harvest in 2008.

Figure 5 Cumulative net biome productivity (NBP) and SOC stock variation modelled by RothC (ΔSOC_RothC) for cover crop (panel a), green dwarf bean (panel b), silage maize (panel c), seed potato (panel d), sugar beet (panel e) and winter wheat (panel f). The cumulative bias between the two variables is shown in the inset. Black circles indicate farmyard manure application. The size of the circles corresponds to the 10th, 50th and 90th percentiles of the fraction of mineralization of the DPM and RPM pools (f_res). n is the number of sowing-to-sowing periods per crop. Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean.

A strong agreement between variations of cumulative NBP and ΔSOC_RothC was found for cover crop, with a cumulative difference of 0.48 ± 0.40 kg C m⁻², which persisted whether farmyard manure was added or not (Fig. 5a). The point at negative cumulative NBP and ΔSOC_RothC variations corresponds to the phacelia – clover crop in 2024, which was truncated in December 2024, resulting in a partial mineralization of C residues after crushing (Fig. 5a). All other cover crop periods had a positive variation of cumulative NBP and ΔSOC_RothC, highlighting SOC stock decrease. Only green dwarf bean showed a cumulative difference of −0.04 ± 0.00 kg C m⁻², but the interpretation is limited with only one crop period in 2021 (Fig. 5b), similarly to maize with a single crop period in 2012 (Fig. 5c).

Winter wheat was an important driver of the discrepancy between cumulative NBP and ΔSOC_RothC over the five-rotation period, with a marked underestimation of C sequestration via NBP compared to RothC over 10 crop periods and a cumulative difference of 1.55 ± 0.54 kg C m⁻² (Fig. 5f). However, low f_res indicated that a significant fraction of the imports of C to the model was not mineralized to conclude that winter wheat significantly sequesters C in the soil. The same observation regarding f_residue applies to all crops except cover crops (Fig. 5b, c, d, e). Since cover crop and winter wheat are part of the same agronomic sequence which aims to cover the soil before the next spring crop (sugar beet or potato), one could group these two crops together to form a unique winter wheat – cover crop sequence spanning from sowing of winter wheat to sowing of the next main crop after cover crop crushing. The sum of the difference between cumulative NBP and ΔSOC_RothC of these new sequences was still driving the positive discrepancy between cumulative NBP and ΔSOC_RothC (0.70 ± 0.63 kg C m⁻², Fig. 6a, b). Importantly, the extension of the sowing-to-sowing period for WW-CC sequences resulted in a higher degree of mineralization (Fig. 6b). Although addition of farmyard manure achieved negative variations of cumulative NBP and ΔSOC_RothC associated with high mineralization rate for two periods in 2009–2010 and 2015–2016, high uncertainties prevent conclusive sink identification for these specific sequences (Fig. 6b). Only one winter wheat – cover crop sequence, corresponding to the phacelia – clover crop in 2024, was a net C sink due to farmyard manure application but still achieved an incomplete mineralization (Fig. 6b). The three remaining winter wheat crop periods alone were the ones in Rotation 1 (2004–2006, 2006–2008 and 2010–2012) when no cover crops were sown after harvest and bare soil was there until sowing of the seed potato, sugar beet and maize respectively (Fig. 6a). For these three winter wheat crops, no significant sink or source pattern, nor bias between cumulative NBP and ΔSOC_RothC could be observed (Fig. 6a).

Figure 6 Cumulative net biome productivity (NBP) and SOC stock variation modelled by RothC (ΔSOC_RothC) for winter wheat without subsequent cover crops (panel a) and for the winter wheat – cover crop sequence (panel b). The cumulative bias between the two variables is shown in the inset. The size of the circles corresponds to the 10th, 50th and 90th percentiles of the fraction of mineralization of the DPM and RPM pools (f_res). Black circles indicate farmyard manure application. n is the number of sowing-to-sowing periods per crop. Abbreviations: WW, winter wheat; CC, cover crop.

4 Discussion

4.1 Discrepancies between RothC simulations and EC

The annual SOC stock loss rate derived from NBP was estimated at approximately 0.07 kg C m⁻² yr⁻¹ on average. At this rate, the initial SOC stock in the 0–60 cm soil profile (5.99 kg C m⁻²) would theoretically be depleted within about 85.6 years, assuming a linear and uninterrupted loss. However, this simple projection should be interpreted with caution, as it does not fully reflect long-term soil carbon dynamics. In particular, the site has been under continuous cultivation for nearly a century, with the current crop rotation maintained for at least 40 years. Under such stable management conditions, SOC dynamics are more likely to approach a quasi-equilibrium state rather than follow a sustained linear decline over time. For instance, Buysse et al. (2013) demonstrated at a long-term experimental site near Lonzée that SOC stock equilibrium is typically achieved after 20 years of constant farmyard manure and crop residue treatments, resulting in minimal subsequent effects on SOC concentration. While these observations suggest that the EC-derived SOC stock loss rate may be overestimated, these results on cumulative NBP lie within the range of previous studies conducted on similar crop rotations using EC data. For instance, Poyda et al. (2019) found an annual loss of 0.12 kg C m⁻² for six intensive crop rotations over a period of eight years in Germany. Emmel et al. (2018) showed that the Oensingen site in Switzerland (site code CH-Oe2) was a net carbon source with a loss of 0.12 kg C m⁻² yr⁻¹ over 13 years for a crop rotation including winter barley, winter wheat, winter rapeseed, peas and potato. Similar results were found across Europe for similar rotations under similar climate (Kutsch et al. (2010): 0.10 kg C m⁻² yr⁻¹, Ceschia et al. (2010): 0.14 kg C m⁻² yr⁻¹, Buysse et al. (2017): 0.09 kg C m⁻² yr⁻¹, Loubet et al. (2011): 0.13 kg C m⁻² yr⁻¹). However, these findings strongly contrast with SOC stock variations derived from regional inventories. For instance, historical soil databases (Goidts and Van Wesemael, 2007) and modelling (Meersmans et al., 2011) have shown that croplands have lost on average 0 to 0.02 kg C m² yr⁻¹ in the plough layer over 50 years in Wallonia, reflecting a limited to negligible decline. Although this trend was supported by RothC simulations and soil inventories at the Lonzée cropland (Fig. 4), our results show large uncertainty bounds in ΔSOC_RothC over the five-rotation period, underscoring the need for cautious interpretation.

Furthermore, the comparison between ΔSOC_RothC and cumulative NBP may be biased by the omission of subsoil carbon dynamics in RothC. While EC fluxes integrate CO₂ exchange across the entire soil profile, RothC simulations only accounted for the top 23 cm of soil (Coleman and Jenkinson, 1996). C inputs from belowground biomass for winter wheat and silage maize were estimated using measured shoot-to-root ratios and harvest index values. These estimates were aligned with the RothC vertical discretization by assigning root-derived inputs to the topsoil layer, under the assumption that most root biomass is concentrated within this horizon (Bolinder et al., 2007). As a result, all carbon inputs were conceptually added into the topsoil and no inputs to the subsoil were considered. Although evidence suggests that most CO₂ efflux in croplands originates from the topsoil (Khan et al., 2023; Wang et al., 2018; Wiaux et al., 2015) and that subsoil SOC stocks are relatively stable due to slow turnover rates (Alcántara et al., 2017) or textural protection (Sierra et al., 2024), recent studies indicate that root–subsoil interactions, particularly in cereals, can substantially contribute to changes in subsoil carbon stocks (Jones et al., 2026; Kätterer et al., 2014; Patra et al., 2022, 2024). In the present study, measurements conducted during pre-harvest field campaign in 2022 showed that winter wheat root biomass in the subsoil (30–100 cm) accounted for 16.5 ± 4.2 % of the total root biomass over the whole soil profile. Therefore, a significant part of C inputs for winter wheat and maize could have been attributed to a subsoil layer, which could have exacerbated the difference between RothC simulations and cumulative NBP by enhancing SOC sequestration due to greater stabilization with depth, thus widening the gap with cumulative NBP. On the other hand, this could also decrease ΔSOC_RothC if priming of older carbon, enhanced mineralization, or lateral exports dominated. To resolve this uncertainty, future work should assess the sensitivity of overall ΔSOC to belowground carbon allocation, depth-dependent turnover rates, and priming effects. Recent modelling frameworks have started to address these processes explicitly – including the implementation of RothC within JULES across a 3 m soil profile (Nakhavali et al., 2018), rCTOOL (Giannini-Kurina et al., 2025), and ORCHIDEE-SOM (Camino-Serrano et al., 2018). While SOC stock measurements for topsoil and subsoil are available at our site and could support such an extension, implementing and validating a multi-layer model was beyond the scope of the present study, which focused on the topsoil carbon budget over observed rotation sequences. Coupling these multi-layer approaches with the long-term eddy covariance and soil inventory dataset presented here therefore constitutes a particularly well-constrained opportunity for future work.

A complementary perspective to further constrain SOC stock dynamics at BE-Lon would be the measurement of functional soil organic matter fractions, such as particulate organic matter (POM) and mineral-associated organic matter (MAOM). Bulk SOC stock measurements integrate pools with contrasting turnover times and stabilization mechanisms, and may therefore mask compensating changes between a labile, plant-derived POM pool and a more persistent, mineral-stabilized MAOM pool (Cotrufo et al., 2019; Lavallee et al., 2020). Partitioning SOC into these measurable fractions would provide a more mechanistic basis for interpreting the near-equilibrium trend observed here and would offer observational constraints that map more directly onto recent process-based models structured around measurable pools (e.g. MEMS, Millennial; Abramoff et al., 2022) than the kinetically defined pools of RothC. Additionally, constraining model uncertainty using supplementary observational data, such as mineralization curves to improve process representation in RothC (Mondini et al., 2017) would support more robust conclusions.

Finally, an important step in the calculation of NBP from NEE and lateral C fluxes was the omission of total dissolved organic and inorganic carbon (DIC-DOC) in Eq. (1), due to a lack of available measurements and limited representation in the model framework. Notably, DOC and DIC fluxes have been shown to be non-negligible in croplands and can represent a significant component of the carbon budget (Nakhavali et al., 2021, 2024). Yearly losses yielded 12.5, 16.0 and 11.2 g C m⁻² yr⁻¹ at the cropland ICOS sites of Grignon, Klingenberg and Carlow, respectively (Kindler et al., 2011). Comparable flux magnitudes have been reported in other croplands (Manninen et al., 2018; Walmsley et al., 2011) as well as in large-scale modelling studies across Europe (Gommet et al., 2022; Janssens et al., 2003). Including DIC and DOC fluxes in the EC-derived C budget would therefore increase the discrepancy between cumulative NBP and ΔSOC_RothC, while leaving unchanged the conclusions of the present study regarding the differences observed among the SOC stock–based methods.

4.2 Uncertainty propagation and its impact on trend significance

The large uncertainty bounds on ΔSOC_RothC indicates that RothC simulations must be interpreted with caution. By the end of the simulation period, the standard deviation of simulated SOC stock variation (0.73 kg C m⁻²) was larger than the corresponding increase in SOC stock (−0.43 kg C m⁻²), underscoring considerable predictive uncertainty on SOC stock predictions. The variability observed here exceeded that reported in previous studies (Filip et al., 2025; Morais et al., 2019; Pesce et al., 2024), likely because, in addition to C inputs and meteorological drivers, we also accounted for a moderate hypothesized uncertainty in the model's functional relationships (20 % on model parameters, Sect. 2.6). However, this study did not consider potential correlations among the perturbed parameters and drivers. All Monte Carlo simulations were based on independent sampling. Several of these variables are however likely correlated: temperature and evapotranspiration typically co-vary within the same climatic events, decomposition rate constants are usually derived from a joint calibration, and biomass inputs may respond to the same climatic drivers used to modify decomposition. Ignoring these dependencies can bias the estimated uncertainty in either direction: positive, reinforcing correlations would lead independent sampling to underestimate true variability, whereas opposing correlations could cause overestimation. As the correlation structure among RothC parameters and drivers was not quantified at our site, we cannot determine a priori whether the reported uncertainty bounds over- or underestimated the true predictive uncertainty. Future applications should consider a joint sampling scheme informed by the empirical covariance among climatic drivers and model parameters. Additionally, using multiple soil C models within an ensemble framework could provide a more robust assessment of uncertainties arising from model structural differences (Bruni et al., 2022; Farina et al., 2021). Recent efforts to develop a single interface for generating projections from multiple SOC models in forest soils pave the way towards more accessible and robust predictions of SOC stocks over croplands (Bruni et al., 2026). Nonetheless, our results suggest that the 2017 SOC stock measured by physical sampling fell within the simulated range, thereby supporting the overall trend of equilibrium of SOC stocks at BE-Lon over the 21-year period. Confirmation of the apparent long-term stability requires additional SOC measurements, ideally at least ten years after the 2017 sampling (Schrumpf et al., 2011), which is required in the ICOS protocols.

Error propagation in EC measurements is critical for reliable carbon budgets, as uncertainties arise from instrumental limitations, missing data and gap-filling. In this study, we quantified uncertainty contributions from key processing steps over the five-rotation period. Random measurement errors contributed minimally to the total error on NEE (0.98 g C m⁻² yr⁻¹), consistent with literature values (0–5 g C m⁻² yr⁻¹) (Du et al., 2014; Vekuri et al., 2023). The u^* threshold uncertainty reached 8.09 g C m⁻² yr⁻¹, which is slightly lower than the range of reported values (10–100 g C m⁻² yr⁻¹) (Papale et al., 2006) likely thanks to sowing-to-harvest seasons in REddyProc that constrained surface roughness variability. Spectral correction uncertainty was comparable at 6.29 g C m⁻² yr⁻¹, reflecting methodological robustness of the spectral correction workflow. These two dominant sources on NEE uncertainty (u^* filtering and spectral correction) together accounted for ∼ 30 % of the total annual uncertainty on cumulative NBP, the other 70 % being shared by the uncertainty on C_imp and C_exp. The total cumulative standard deviation on NEE represented a very small fraction of the cumulative flux (∼ 3 % on the −7.72 kg C m⁻²), thereby indicating that the post-processing workflow of EC data minimized the uncertainty associated with flux estimation and ensured a high level of confidence in the derived cumulative carbon balance.

4.3 Potential underestimation of C sequestration

This study showed a systematic overestimation of cumulative NBP when compared to ΔSOC_RothC, which was observed for winter wheat, maize, seed potato and sugar beet (Fig. 5), and suggested a potential underestimation of C sequestration by the EC technique. To further investigate the mismatch between the flux-based approach and the SOC modelling approach, we calculated the carbon balance ratio (CBR), which is defined as the ratio of the integrated NEE to the difference between the carbon exported and imported plus the difference between final and initial SOC (Baker and Griffis, 2005). In the absence of SOC stock changes, the vertical CO₂ flux balance should match the net lateral carbon balance, leading to a CBR of ∼ 1. At the Lonzée site, the CBR over five full rotations was ∼ 82.5 %, suggesting that EC measurements likely captured only a part of the true CO₂ exchange over the cropland, assuming no systematic biases affected lateral carbon fluxes or inventory estimates. This finding for CO₂ fluxes is consistent with the well-documented non-closure of the surface energy balance, where the sum of turbulent sensible and latent heat fluxes is likely smaller than the available energy (Aubinet et al., 2012; Foken, 2008; Oncley et al., 2007; Twine et al., 2000). The energy balance ratio (EBR) at BE-Lon was 60.8 % over the five-rotation period (Arquin, pers. comm.), which falls within the lower range of values typically reported for croplands (Nicolini et al., 2026; Stoy et al., 2013). One commonly proposed explanation is the contribution of large-scale turbulent structures with timescales exceeding the EC averaging period, which generate dispersive fluxes that are not captured by single-tower measurements (Foken et al., 2011; Gao et al., 2019; Mauder et al., 2020; Wohlfahrt and Widmoser, 2013). It has also been suggested that missing terms in the classical energy balance equation may contribute to the observed imbalance (Durand, 2022; Sun et al., 2021; Sun and Tribbia, 2025). To test whether the discrepancy between the EC-based C budget, the model estimates and the soil C inventories could arise from a fraction of the turbulent CO₂ exchange being missed by EC at the ecosystem surface, we computed the EBR following Volk et al. (2021) and attributed a fraction (α) of the energy balance non-closure to NEE, varying α between 0 % and 100 %. The optimal α was identified by comparing the resulting gain in NEE after correction (ΔNEE) with cumulative NBP. α=52 % yielded a cumulative ΔNEE of the same magnitude as cumulative NBP (1.45 kg C m⁻²), indicating that correcting NEE for only part of the energy balance non-closure may help reconcile EC, RothC and soil C inventories. This supports the hypothesis that eddy covariance turbulent fluxes are underestimated in absolute value at BE-Lon, while also implying that the degree of underestimation differs between turbulent energy fluxes and CO₂ fluxes. Whether energy balance non-closure can be systematically translated into NEE corrections nevertheless remains uncertain (Mauder et al., 2024; Zhou et al., 2024). These findings therefore call for caution when interpreting absolute NEE-derived carbon budgets, while highlighting the value of independent model-based approaches, such as RothC, to benchmark and contextualize EC measurements.

4.4 The role of amendments and cover crops on SOC sequestration

This study demonstrates that farmyard manure application during a winter wheat – cover crop sequence led to a significant increase in SOC stocks, as consistently indicated by both cumulative NBP and ΔSOC_RothC variations under high residue mineralization conditions (Fig. 6b). In contrast, without manure amendment, the additional carbon inputs from cover crop incorporation and winter wheat harvest residues appeared insufficient to offset heterotrophic respiration losses, resulting in a moderate net carbon loss (Fig. 6b). Manure applications can enhance SOC stocks through two non-exclusive pathways: (i) by boosting crop productivity and cover crop growth, thereby increasing the amount of plant residues and rhizodeposition returned to the soil (Stewart et al., 2007); and (ii) by directly supplying organic matter with a composition and nutrient content that promote the stabilization of a larger fraction of input C in soil organic matter compared to crop residues alone, as suggested by global syntheses of manure effects on SOC (Maillard and Angers, 2014). These findings highlight the critical role of organic amendments in enhancing the carbon sequestration potential of the BE-Lon cropland and raise important questions regarding the capacity of cover crops alone – in the absence of external organic inputs – to drive net SOC accumulation under the prevailing climatic and management conditions.

The potential to further enhance carbon sequestration through increased manure application depends on factors such as soil and climate adaptability, economic feasibility for farmers, impacts on food production, and the extent to which these practices are already implemented. It is well known that manure applications have been shown to increase SOC stocks (Bolinder et al., 2020; Poulton et al., 2018; Sail et al., 2025), making organic amendments a promising strategy to increase SOC stocks in agricultural systems (Bai et al., 2023; Chenu et al., 2019). However, the potential for further SOC gains from increasing manure application depends on the soil's capacity to stabilize added carbon. This capacity is determined by the number of mineral binding sites, which is closely related to the proportion of fine silt and clay particles (Carter et al., 2003). Soils are therefore expected to have a finite carbon storage capacity in aggregate fractions, with an estimated upper limit of approximately 50 g C kg⁻¹ for European soils (Cotrufo et al., 2019). Consequently, soils with low initial SOC content in the fine fraction typically have high carbon sequestration potential due to a greater capacity for long-term stabilization in association with the mineral phase (Georgiou et al., 2022; Gross and Glaser, 2021). However, recent studies indicate that soils with high initial SOC can continue to sequester carbon when carbon inputs increase (Orgill et al., 2017). It could be expected that increasing manure applications leads to measurable gains in SOC stocks (Joona et al., 2024; Koga, 2017; Liang et al., 2023; Shakoor et al., 2021).

The role of cover crops in enhancing the carbon sink capacity of croplands remains debated in the literature, with studies reporting a wide range of outcomes from net carbon gains to neutral or slightly negative balances. For instance, a long-term experiment in France demonstrated that residue incorporation from cover crops was largely compensated by the associated CO₂ respiration flux following crushing, in a cropland with a SOC content comparable to BE-Lon (Breil et al., 2023). Meta-analyses notably consistently report SOC stock gains under cover crops relative to bare soils (McClelland et al., 2021; Poeplau and Don, 2015; Scavo et al., 2022). Process-based model simulations in Flanders (Spotorno et al., 2024) or Germany (Seitz et al., 2023) similarly suggest beneficial effects on SOC stocks. However, other studies point to more limited or inconsistent outcomes. Blanco-Canqui (2022) reported only minor effects of cover crops on SOC sequestration, a 65-year long-term experiment in Belgium found marginal SOC increases attributable to cover crops (Sail et al., 2025), and Schepens et al. (2025) concluded that winter cover crops do not systematically increase SOC stocks across Northwestern Europe. The extent to which cover crops enhance the carbon sink capacity of croplands is therefore site-specific, and depends on factors including the duration of cover crop integration, initial SOC content, soil texture, climate and, critically, plant species (Blanco-Canqui, 2022; Breil et al., 2023; Zhang et al., 2022). At BE-Lon, mustard accounted for more than 80 % of cover crop periods, with oat and faba bean mixtures comprising the remainder. Mustard appears to have limited SOC sequestration potential (Antosh et al., 2020), in contrast to ryegrass – also used in Western Europe (Smit et al., 2019) – which can substantially enhance SOC stocks (Jensen et al., 2022; Poeplau et al., 2015). The dominance of a low-sequestration cover crop species at the BE-Lon site may therefore partly explain the absence of a net carbon sink effect from cover cropping. Furthermore, incorporating N-rich cover crops like faba beans may lead to increased N₂O emissions offsetting the carbon storage benefits (Li et al., 2024).

This study could not confirm whether winter wheat or maize as contributors to SOC sequestration, owing to the incomplete mineralization of the DPM and RPM pools following harvest (Fig. 5). When winter wheat data were pooled with the subsequent cover crop periods to achieve higher DPM mineralization rates over a unique agronomic sequence, both cumulative NBP and ΔSOC_RothC consistently indicated a moderate net carbon source trend – except when farmyard manure was applied (Fig. 6). This pattern is consistent with the broader finding that organic amendments, rather than crop rotation design alone, appear to be the primary lever for SOC accumulation at BE-Lon. While the inclusion of cereal crops with deep root systems has been proposed as a strategy to increase belowground C inputs and promote SOC gains in many studies (Patra et al., 2022, 2024; Poeplau et al., 2021; Sun et al., 2018), this effect may be insufficient to counteract heterotrophic respiration losses under the management and climatic conditions prevailing at this site. Under the conditions of this study, modifying the rotation to increase the proportion of cereal crops does not appear to be a sufficient strategy for improving net SOC sequestration.

Although this study showed that manure application was a key driver of SOC accumulation, the associated increase in greenhouse gas emissions, particularly nitrous oxide (N₂O), may partially or even fully offset the climate benefits of enhanced carbon sequestration (Guenet et al., 2021). In general, manure amendments lead to significantly higher N₂O emissions than mineral nitrogen fertilization alone, with the magnitude of this effect strongly influenced by the $\mathrm{C}/\mathrm{N}$ ratio of the applied material (Zhou et al., 2017). While the present study did not quantify N₂O emissions following manure application, additional inputs, even within the limits set by the European Union Nitrates Directive (91/676/EEC; European Commission, Directorate-General for Environment, 2000) and farm-level agronomic constraints, could exacerbate ammonia volatilization, nitrate leaching, and N₂O emissions, thereby reducing or negating the net greenhouse gas mitigation potential at the cropland scale (Maenhout et al., 2024). Moreover, manure use is constrained by broader environmental concerns, including water contamination, nitrate leaching, the spread of antibiotic resistance, and risks to ecosystems functioning (Zhen et al., 2020, Wang et al., 2019). Overall, the net climate mitigation potential of manure application is likely limited, as many of its benefits reflect redistribution of existing carbon stocks rather than direct sequestration of atmospheric carbon into the soil (Powlson et al., 2011, 2012).

5 Conclusion

This study compared decadal soil inventories, RothC simulations, and EC measurements with lateral C flux inventories to assess SOC stock changes at BE-Lon across five crop rotations. EC-derived carbon budgets indicated high SOC losses inconsistent with the long-term agricultural and management history of the site, while RothC simulations and soil inventories supported near-equilibrium dynamics. Uncertainties associated with spectral correction and friction velocity threshold determination were comparable in magnitude to the uncertainty in lateral C flux estimates within the EC-based budget. However, all uncertainties remained smaller than the cumulative NBP signal, enabling the identification of a strong source trend. Comparison between cumulative NBP and RothC simulations of SOC stock variation revealed a systematic positive bias for most crops, likely reflecting an underestimation of C sequestration by the EC technique, which is consistent with the well-documented energy balance non-closure at cropland sites. Both approaches identified the winter wheat – cover crop sequence as a moderate net C source in the absence of organic amendments. Farmyard manure application emerged as the primary driver of SOC accumulation at the site, as shown by the convergence between the two approaches under high DPM mineralization. However, evaluating the true mitigation potential of this practice requires accounting for potential trade-offs with N₂O emissions in the full greenhouse gas balance. More broadly, this study highlights the value of combining process-based models, multi-decadal soil inventories, and eddy covariance flux tower to constrain SOC budgets, and improve understanding of SOC stock drivers on croplands.

Appendix A

Figure A1 Spatial design of soil sampling in 2007 and 2017. The black points indicate 2007 samples location. The blue diamonds indicate 2017 samples locations labelled as “SP-I”. Strata were determined based on Loubet et al. (2026).

Table A1 Soil characteristics at BE-Lon across the vertical profile: carbonates (CaCO₃), soil organic carbon content (SOC), organic matter (OM), nitrogen content (N), bulk density (BD), pH and granulometry. Samples were collected during the field campaign in 2017 following ICOS protocols.

Table A2 Seed density (SD), dry matter (DM_H), carbon content ([C]_H) of the harvested product and the corresponding export of carbon at harvest (C_exp) for the main crop periods. For abbreviations, WW: winter wheat, SP: seed potato, SB: sugar beet, MA: silage maize, GDB: green dwarf bean. All uncertainties are standard deviation.

Table A3 Cumulative values of net ecosystem exchange (NEE), imports of carbon (C_imp), exports of carbon at harvest (C_exp), net biome productivity (NBP) and aboveground crop residues at harvest (CR). Crops are mentioned with the corresponding year. All uncertainties are cumulative standard deviation. Bold text for C_imp highlights the application of farmyard manure. Crop-periods are defined as from sowing to sowing of the next crop. Abbreviations: WW, winter wheat; SP, seed potato; SB, sugar beet; CC, cover crop; MA, silage maize; GDB, green dwarf bean.