Regional assessment and uncertainty analysis of carbon and nitrogen balances at cropland scale using the ecosystem model LandscapeDNDC

The assessment of cropland carbon and nitrogen (C and N) balances plays a key role in identifying cost-effective mitigation measures to combat climate change and reduce environmental pollution. In this paper, a biogeochemical modelling approach is adopted to assess all C and N fluxes in a regional cropland ecosystem of Thessaly, Greece. Additionally, the estimation and quantification of the modelling uncertainty in the regional inventory are realized through the propagation of parameter distributions through the model, leading to result distributions for modelling estimations. The model was applied to a regional dataset of approximately 1000 polygons, deploying model initializations and crop rotations for the five major crop cultivations and for a time span of 8 years. The full statistical analysis on modelling results (including the uncertainty ranges given as ± values) yields for the C balance carbon input fluxes into the soil of 12.4 ± 1.4 t C ha⁻¹ yr⁻¹ and output fluxes of 11.9 ± 1.3 t C ha⁻¹ yr⁻¹, with a resulting average carbon sequestration of 0.5 ± 0.3 t C ha⁻¹ yr⁻¹. The averaged N influx was 212.3 ± 9.1 kg N ha⁻¹ yr⁻¹, while outfluxes of 198.3 ± 11.2 kg N ha⁻¹ yr⁻¹ were estimated on average. The net N accumulation into the soil nitrogen pools was estimated to be 14.0 ± 2.1 kg N ha⁻¹ yr⁻¹. The N outflux consists of gaseous N fluxes composed of N₂O emissions of 2.6 ± 0.8 kg N₂O–N ha⁻¹ yr⁻¹, NO emissions of 3.2 ± 1.5 kg NO–N ha⁻¹ yr⁻¹, N₂ emissions of 15.5 ± 7.0 kg N₂–N ha⁻¹ yr⁻¹ and NH₃ emissions of 34.0 ± 6.7 kg NH₃–N ha⁻¹ yr⁻¹, as well as aquatic N fluxes (only nitrate leaching into surface waters) of 14.1 ± 4.5 kg NO₃–N ha⁻¹ yr⁻¹ and N fluxes of N removed from the fields in yields, straw and feed of 128.8 ± 8.5 kg N ha⁻¹ yr⁻¹.

1 Introduction

Food security and agricultural productivity depend, to a major extend, on the applied nitrogen (N) fertilizers (Klatt et al., 2015). Worldwide, the N fertilizer use for the years 1960 to 2005 increased from 30×10⁶ to 154×10⁶ t (IFADATA, 2015). In Europe, the increase in yields of arable land and grassland systems was 45 %–70 % from 1950 (European Fertilizer Manufacturers Association, 2009) due to the intensification of agricultural production systems. Excessive use of N fertilizers, though beneficially affecting the yield, could have a harmful impact on the environment, e.g. increased gaseous emissions and aquatic fluxes of nitrous oxide (N₂O) to the atmosphere and leaching of nitrate (NO₃) into waterbodies (Erisman et al., 2011; Galloway et al., 2013; Kim et al., 2015)

N₂O poses a twofold environmental threat. On the one hand, it is a strong greenhouse gas with a warming potential that is 300 times greater (in a 100-year time period) than carbon dioxide (CO₂), and on the other hand, it is a major driver of ozone depletion in the stratosphere (Ravishankara et al., 2009). Fertilizer use aimed at an increase in agricultural production is the most crucial anthropogenic source of atmospheric N₂O, which, at present, contributes approximately 45 % of the total anthropogenic N₂O emissions on a global scale (Jones et al., 2014). Because of global population growth and, thus, growing food demand (Godfray et al., 2010), fertilizer use will probably increase. Consequently, the predictions of the current business-as-usual scenarios show doubled anthropogenic N₂O emissions by the year 2050 (Davidson and Kanter, 2014). The European countries have recently set up bilateral agreements in order to reduce N₂O emissions from cultivated croplands (EU-Commission, 2014). Similarly, the European Nitrates Directive (EU-Commission, 2019; Musacchio et al., 2020) aims to reduce NO₃ leaching to waterbodies to avoid an increase in both eutrophication (Camargo and Alonso, 2006) and drinking-water pollution. Because of the hazardous N₂O and NO₃ effects, it is necessary to evaluate agricultural systems in terms of their profitability and productivity and in terms of their impacts on the environment.

The N₂O and NO₃ production and consumption in agricultural lands are regulated, to a large extent, by N plant uptake and, also, by the microbial processes of denitrification and nitrification (Butterbach-Bahl et al., 2013). The factors controlling both the microbial metabolism and plant N uptake are (a) soil conditions (Butterbach-Bahl et al., 2013) and (b) cultivation management practices, e.g. crop rotation, fertilizing amount and timing, and ploughing (Smith et al., 2008).

In order to reach a minimization of the environmental footprint of agricultural production while securing global food security (Garnett et al., 2013), it is mandatory to tighten the N cycling of intensified agricultural systems, e.g. by harmonizing the N demand of crops with soil N availability by means of N fertilization.

Full nitrogen balance inventories provide a comprehensive description of the different N input and output fluxes within an arable system to the scientific community, farmers and policymakers. The assessment of the N balance is essential to optimize nitrogen use and production and to minimize environmental impacts and pollution. In particular, policymaking and regulatory bodies require accurate and robust information on all different nitrogen fluxes to develop effective strategies in agricultural N management. Until now, our descriptions of N cycling in arable land have lacked observations of the full N balance as only a few studies have tried to quantify the total N balance of agricultural systems, e.g. the study by Zistl-Schlingmann et al. (2020) using stable isotope techniques or the study by Schroeck et al. (2019) using process-based modelling.

A recent opinion paper by a large group of leading scientists in the field of process-based ecosystem modelling, Grosz et al. (2023), identified the lack of knowledge of the full N balance and the fact “the scarcity of complete modeled N balances in the literature stems from the reluctance of the scientific community to support the publication of unvalidated modeled results, especially given that the simulation results of these neglected N pools and fluxes may be unrealistic”. It is further noted that “This self-censorship of authors has resulted in a missed opportunity to share knowledge and improve our understanding of modelled processes”.

Grosz et al. (2023) conclude that “including the entire N balance and related parameters should become standard when publishing the results of N model studies”. Grosz et al. (2023) emphasize that this would allow us to assess the robustness of modelled N fluxes and full N balances and to illustrate the diversity and uncertainty of the different process-based modelling approaches, e.g. modelling denitrification processes in soils.

In this analysis, the process-based bio-geochemical model LandscapeDNDC (Haas et al., 2013) was applied to the agricultural cropland systems in the region of Thessaly (Greece). The objective of our study was threefold.

We assess and report the cropland C and N balances, including all associated fluxes such as CO₂, N₂O and NH₃ emissions; NO₃ leaching; and the soil carbon stock changes as demanded by Grosz et al. (2023).

We increase the robustness and trustworthiness of the balance modelling by assessing and quantifying the modelling uncertainty of the simulated C and N balances and flux estimations, as requested before by the IPCC (IPCC, 2019).

We present a regional uncertainty assessment methodology for C and N cycling to advance the balance modelling by propagating 500 joint parameter and input data distributions through the model (each representing a full regional C and N balance inventory simulation), yielding regional result distributions for any modelling estimations.

2 Material and methods

2.1 Model description

LandscapeDNDC is a modular, process-based ecosystem model for simulating the bio-geochemical change in C and N in cropland, forest and grassland systems at both the site scale and regional scale. The modules combined are about plant growth, micro-meteorology, water cycling, physico-chemical plant and microbial C and N cycling, and the exchange processes of the atmosphere and hydrosphere of terrestrial ecosystems. LandscapeDNDC is a generality of the plant development and soil biogeochemistry of the agricultural DNDC and Forest-DNDC (Li, 2000). There has been a successful application of earlier model versions in a number of studies, e.g. for water balance (Grote et al., 2009; Holst et al., 2010), for plant growth (Cameron et al., 2013; Werner et al., 2012), for NO₃ leaching (Thomas et al., 2016; Kim et al., 2015), and for soil respiration and gas emission traces (Chirinda et al., 2011; Kraus et al., 2014; Molina-Herrera et al., 2015). For the initialization of LandscapeDNDC, physical and chemical site-specific soil profile information is used (specified for different soil depths): soil organic carbon (SOC) and nitrogen (SON) content, soil texture (clay, sand and silt content), the plant growth and soil biogeochemistry, bulk density, pH value, saturated hydraulic conductivity, field capacity, and wilting point. Daily or hourly climate data of air temperature (max, min and average), N deposition, precipitation and atmospheric CO₂ concentration are used in LandscapeDNDC in combination with agricultural management practices; e.g. crop planting and harvesting, fertilization (synthetic and organic), or feed cutting and tilling are used to drive LandscapeDNDC simulations. Regarding fertilization management, three types of mineral fertilizers, namely urea, compound fertilizers based on NH₄ and NO₃, and organic amendments (i.e. green manure, farmyard manure, slurry, straw, bean cake and compost), are currently considered. The growth of crops and grasses is similar to the DNDC approach using two major parameters that describe seasonal plant development (cumulative-temperature degree days) and maximum reachable biomass under optimum conditions (Li, 2000), with daily growth limitations due to water and nutrient availability considered. Model parameters describing soil and vegetation characteristics are obtained from an external parameter library. In LandscapeDNDC, the parameterization of the main cultivated commodity crops in Europe is performed using default parameter sets representing an average plant type, while process parameter values for micro-meteorology, water cycle and bio-geochemical processes were obtained from previous validation studies (e.g. Molina-Herrera et al., 2016; Klatt et al., 2015; Rahn et al., 2012), proving that the LandscapeDNDC model could be universally applicable to similar conditions.

For all simulations in the current study, site-specific crop parameterizations were derived in a preceding analysis of various site-scale simulations and validations of yield characteristics across the region. An overview of the crops cultivated at the different study sites and detailed information on the specific crop rotations used to simulate crop growth are provided in Table S2 in the Supplement.

2.2 Case study description and input data

The region of Thessaly is located in central Greece, covering a total area of 14 000 km², where 5000 km² is lowland and approximately 2300 and 6500 km² are semi-mountainous and mountainous land respectively. The plain of Thessaly is considered to be among the largest agricultural lands of the country (Kalivas et al., 2001), accounting for almost 410 000 ha, of which about 370 000 ha is arable land, where almost 80 % is covered by annual crops, and 10 % is covered by perennial crops (ELSTAT: https://www.statistics.gr/el/statistics/-/publication/SPG06/2012, last access: March 2017, 2021.). The crop or plant production of the region is around 14.2 % (ELSTAT, 2012) of the total production of the country (second in Greece).

Soil input data for the region were available from the European project NitroEurope IP (Sutton et al., 2013) based on the European Soil Database (ESDB, 2004) containing soil type and soil profile descriptions of bulk density, SOC content, texture (sand, silt clay), pH value, stone fraction, saturated hydraulic conductivity, wilting point and water-holding capacity for various soil strata (Cameron et al., 2013). A regional soil dataset for the area of interest contained about 1500 spatial polygons, out of which approximately 1000 covered the cultivated cropland that was finally simulated. The climate data for the regional simulations were derived at polygon level from gridded ERA5 climate data for Greece.

2.3 Agricultural management and model input data processing

The total cultivated area and the yields for the years 2010 to 2016, as used in the current analysis, were obtained from the Hellenic Statistical Authority (ELSTAT). Moreover, data associated with the animal capital for the indicated years were also provided (ELSTAT) in order to estimate the annual manure production in the region; however, no data are available on whether the manure is used in croplands and, if this is the case, how much of it is used. For the water management, the percentage of irrigated and non-irrigated land (estimated to be almost 50 % for each case) was also given (ELSTAT), while indicative sets of irrigation management data were acquired through the River Basin Management Plans of the Special Secretariat for Water – Ministry of Environment, Energy and Climate Change (YPEKA, Portmann et al., 2010). The irrigation water volumes were estimated based on the crop needs and the minimum and maximum quantities necessary (according to the literature) while using upscaling tools to obtain the regional values. The fertilization datasets were provided by the Fertilizer Producers and Merchandiser Association (FPMA) for recent years (2010–2016) and are equated to the national-scale annually consumed quantities scaled down to a regional level based on crop patterns in Thessaly's cultivated land.

In this study, the five main crops of maize, wheat, clover, cotton and barley were considered, covering the majority of the cultivated arable land in the region (over 95 %), while the remaining cropland was included, acquiring the final corrected land–crop coverage. In Table 1, the resulting crop rotation scenarios (R1–R5) are presented for the evaluation period of 2012–2016. Note that each rotation sequence (R1–R5) is shifted in time such that, for each year, each crop appears in one rotation. Based on the crop cover contribution in each simulated year, the crop rotation contribution factors were estimated and are summarized in Table 2. The management practices were based on the general agricultural practices applied in the region and on information provided by farmers.

Table 1 Summary of the crop rotation scenarios (R1–R5) for the region of Thessaly. The crop abbreviations corn, wiwh, clover, cott and wbar refer to maize (food corn and silage maize), winter wheat, clover (legume feed crops such as alfalfa or vetch), cotton and winter barley respectively.

Table 2 Crop cultivation area contribution per year to the aggregation of the five rotations; data constant across the region of Thessaly.

2.4 Uncertainty analysis

As stated in the IPCC 2006 guidelines, updated in 2019, the assessment of uncertainty is considered to be a major and crucial or mandatory component when compiling regional or national GHG emission inventories (Larocque et al., 2008). The difference in terms of the scale at which the model is used results in divergent errors in the C and N dynamic predictions across different climate zones and scales. Thus, uncertainty analysis is a crucial step towards a higher-quality decision-making process. The sources of uncertainty can vary and are related to (a) the initial conditions (starting values), (b) the drivers (e.g. climate and crop management data), (c) the conceptual model uncertainty and (d) the parameter uncertainty of the various processes (Refsgaard et al., 2007; Wang and Chen, 2012).

Santabárbara (2019) performed a Bayesian model calibration and uncertainty analysis using a Monte Carlo Markov chain (MCMC) approach targeting uncertainties associated with the data (bulk density, SOC, pH, clay content) of the initial soil conditions, drivers (cropland management such as fertilization or manure rates and timing, harvest and seeding timing, tillage timing) and bio-geochemical process parameterizations.

In order to identify the most sensitive process parameters with a reduced number of model simulations, the Morris method (Morris, 1991) obtains a hierarchy of parameter influences on a given output (gaseous N fluxes) and evaluates whether or not a non-linearity exists. Morris (1991) proposed that this order can be assessed through the statistical analysis of the changes in the model output that are produced by the one-step-at-a-time changes in n number of proposed parameters. Incremental steps in the number of parameters lead to the identification of parameters having substantial influence over the concerned results without neglecting that some effects could cancel each other out (Saltelli et al., 2000), leading to the identification of the 24 most sensitive process parameters (Myrgiotis et al., 2018b; Houska et al., 2017).

2.5 Metropolis – Hastings algorithm

The Markov Chain Monte carlo (MCMC) Metropolis–Hastings algorithm results in numerous parameter sets that approximate the posterior joint parameter distribution by performing a random walk through the space of joint parameter values. This probability evaluation of the data obtained from each step leads to the update of the initial uniform parameter distributions. Bayes' formula relating conditional probabilities may become a powerful and practical computational tool when combined with Markov chain processes and Monte Carlo methods or the so-called Markov Chain Monte carlo (MCMC) approach. A Markov chain is a special type of discrete stochastic process wherein the probability of an event depends only on the event that immediately precedes it. Integrating parameters (θ) and observation data (D) into Bayes' rule results in the following formula:

$$\begin{array}{}\text{(1)}& P\left(\mathit{\theta }|D\right)={\displaystyle \frac{P\left(D|\mathit{\theta }\right)\cdot P\left(\mathit{\theta }\right)}{P\left(D\right)}},\end{array}$$

where P(D|θ), the probability of the data, is used to obtain the probability of these parameters updated by the data, P(θ|D), where the evidence is computed as

$$\begin{array}{}\text{(2)}& P\left(D\right)=\int \text{likelihood}\cdot \text{prior}\cdot d\mathit{\theta },\end{array}$$

where P(D) can be numerically approximated with the aforementioned MCMC method (Robert and Casella, 2011).

The method uses prior knowledge concerning the sources of the model uncertainty to obtain a narrowed posterior distribution for each one of the sources. By propagating the parameter distributions through the model, the overall uncertainty in the model results can be quantified.

In a previous study by Santabárbara (2019), an extensive sensitivity analysis on all soil bio-geochemical process parameters, soil initial data and arable management data was performed, identifying the 24 most sensitive process parameters (listed in the Supplement), the most sensitive soil initial data (soil profile data on bulk density, soil organic carbon content and pH value) and the most sensitive management information (fertilization and manure N rates, tilling depth) in relation to aquatic and gaseous N fluxes from arable soils. This was depicted in the MCMC simulation through sampling a combination of 24 parameter values, 3 values of soil initial data and 3 items of management information. The sampling of the soil initial data and the management data was performed as perturbations to the existing data. For each quantity, a perturbation was sampled individually and applied to all corresponding values in the soil profile or to all years in the management description. The MCMC simulation performed by Santabárbara (2019) simulated more than 100 000 iterations for various arable sites until the MCMC simulation converged towards a stable combined posterior distribution of parameter values and soil and management input data perturbations. In the current analysis, we have sampled 500 joint parameter and input data perturbation sets from the posterior distributions reported by Santabárbara (2019), and we deployed them in simulations (propagation through the model) for the regional inventory, leading to 500 inventory simulations. Afterwards, a statistical analysis was applied to estimate the updated regional and temporal result distributions.

2.6 Statistical methods and data aggregation

2.6.1 Regional result aggregation

One full regional inventory simulation consists of 10 individual inventory simulations: five different crop rotations for irrigated conditions and five different crop rotations for rain-fed conditions were simulated in parallel (see Sect. 2.3). The results of the crop rotations were aggregated according to the crop shares per year (see Table 2), accounting for all the effects of the different crops cultivated in the region under irrigated and rain-fed conditions. The final inventory simulation results were obtained by considering irrigated versus rain-fed water management. The final inventory contains simulation results aggregated to area-weighted yearly means across the total simulation domain, accounting for the cropland area of each polygon.

2.6.2 Uncertainty quantification and statistical analysis

A regional aggregation was performed for all 500 uncertainty simulations. All the uncertainty results were finally reported via statistical measures evaluating the 500 regional uncertainty simulation runs reporting mean values, standard deviations, medians, and the 25th and 75th interquartile ranges (IQRs, Q25 to Q75).

3 Results analysis and evaluation

The simulation time span was from 2009 to 2016, while the years 2009–2011 were used as spin-up to get all soil C and N pools into equilibrium after the initialization. Therefore, reported simulation results are limited to the years 2012–2016. The assessment of the regional C and N balances (CB and NB) were obtained and, as a consequence of the uncertainty quantification, resulted in distributions; they were therefore reported by statistical measures such as the mean and/or median or interquartile ranges of the uncertainty ensemble.

3.1 Regional yield simulations and validation

The evaluation of the model performance in estimating the NB and CB components was analysed based on the comparison of the simulated yield values with the observed yield data provided by the Hellenic Statistical Authority (ELSTAT), averaged for the total simulated period.

Crop yields and feed production

For model validation, datasets of crop yields from the Hellenic Statistical Authority (ELSTAT) were used. Table 3 summarizes the aggregated regional crop yields for all the simulated years and the respective mean, median and standard deviation values that resulted from the statistical analysis of the simulation results together with the observed yield and feed production provided by the Hellenic Statistical Authority (ELSTAT). Simulated yields consist of the cotton bolls of cotton; the clover feed, which is the total cut and harvested aboveground biomass; the grain yield of wheat and barley; and the grain ear and the stems of maize. Based on the observations, maize appears to be the dominant crop, with an average yield of 12 t ha⁻¹, followed by clover products at 8.4 t ha⁻¹. The rest of the three crop yields appear to be of the same order of magnitude from 3.3 up to 3.4 t ha⁻¹.

Table 3 Simulated and observed yields and feed production (t dry matter ha⁻¹) in the region of Thessaly. All results are based on statistical aggregation across all polygons, rotations and years and, finally, across all 500 UA (uncertainty analysis) inventory simulations. The observed values of dry matter (DM) are provided by the Hellenic Statistical Authority.

^* Observation data for maize did not distinguish between food corn and silage maize.

Additionally, the simulated average yield of cotton was estimated to be 3.3 ± 0.8 t DM ha⁻¹, the simulated average yield of wheat was estimated to be 3.6 ± 0.9 t DM ha⁻¹, the simulated average yield of barley was estimated to be 4.5 ± 1 t DM ha⁻¹, and the simulated average yield of maize was estimated to be 9.9 ± 1.4 t DM ha⁻¹. As for the feed, the clover was estimated to be 9.6 ± 0.6 t DM ha⁻¹. The average nitrogen use efficiency (NUE) across time and space is 63.29 %.

Figure 1 presents the uncertainties of the simulated crop yield across the whole evaluation time span of 2012–2016 in both irrigated and rain-fed conditions. As shown, corn shows a much more narrow distribution with a higher median for the irrigated scenario compared to the rain feed. On the contrary, winter barley has a wider distribution and a slightly higher median for the irrigated scenario and, also, a wider extreme value variation. As for cotton, the distribution appears to be bimodal for the rain-fed scenario, in which the median is also lower than the one in the irrigated scenario. In addition, the extreme value variation is wider in the latter case. Finally, for the example of winter wheat, the irrigated and rain-fed scenarios reach the same results.

Figure 1 Simulated crop yield uncertainties across the evaluation time span of 2012–2016 for (a) irrigated and (b) rain-fed conditions. Horizontal lines indicate the median, mean, maximum and minimum values of the distributions.

3.2 Regional carbon and nitrogen balance

3.2.1 Carbon balance (CB)

For the CB, Fig. 2 presents average C input fluxes into the soil of 12.4 ± 1.4 t C ha⁻¹ yr⁻¹ (with an interquartile range (IQR) (from Q25 to Q75) of 12.1 to 13.2 t C ha⁻¹ yr⁻¹) and output fluxes of 11.9 ± 1.3 t C ha⁻¹ yr⁻¹ (with an IQR from 11.6 to 12.7 t C ha⁻¹ yr⁻¹). The resulting carbon sequestration was estimated to be 0.5 ± 0.3 t C ha⁻¹ yr⁻¹, with an IQR from 0.4 to 0.7 t C ha⁻¹ yr⁻¹ (data summarized in Table 4).

Figure 2 Carbon balance for cropland cultivation for the region of Thessaly: (a) averaged carbon balance (in t C ha⁻¹), (b) total carbon balance of cropland soils (in mio. t C) and (c) averaged fluxes across the region and the years 2012–2016. Positive change indicates soil C sequestration.

The input fluxes consist of annual gross primary productivity (GPP) of 11.7 ± 1.4 t C ha⁻¹ yr⁻¹, with an IQR from 11.4 to 12.4 t C ha⁻¹ yr⁻¹, and the carbon in manure applied to soils is estimated to be 0.7 ± 0.001 t C ha⁻¹ yr⁻¹ (see Table 4). This compares, on the other hand, to respirative carbon fluxes from the soil to the atmosphere (TER) of 8.5 ± 1.1 t C ha⁻¹ yr⁻¹, with an IQR from 8.2 to 9.1 t C ha⁻¹ yr⁻¹, and carbon fluxes via exported crop yields and feed (including all straws and removed crop residues) of 3.4 ± 0.3 t C ha⁻¹ yr⁻¹, with an IQR from 3.4 to 3.6 t C ha⁻¹ yr⁻¹. The aggregation of the carbon fluxes to the regional level of approximately 360 000 ha of cropland results in 4.25 ± 0.20Ṁ t C yr⁻¹ in terms of GPP and 0.25 ± 0.01 M t C yr⁻¹ of carbon influx via organic fertilizers compared to 3.08 ± 2.97 M t C yr⁻¹ of TER and 1.24 ± 0.05 M t C yr⁻¹ of carbon exports via crop yields and feed production, leading to a net carbon sequestration of 0.5 ± 0.3 M t C ha⁻¹ yr⁻¹, with an IQR from 0.4 to 0.7 M t C ha⁻¹ yr⁻¹ (M t C indicates million tonnes of carbon).

Table 4 Carbon balance (per hectare) assessment and uncertainty analysis of the of cropland cultivation in the region of Thessaly, Greece. The ¹ mean; ² standard deviation; ³ median; ⁴ Q25, 25th quartile; and ⁵ Q75, 75th quartile are applied across the 500 values for the quantities in this table. ⁶ C inputs as the sum of the absolute values of all the input fluxes of the 500 simulations. ⁷ C outputs as the sum of the absolute values of all the output fluxes of the 500 simulations. ⁸ SOC changes as the difference between the input and output fluxes of each of the 500 simulations. Note that the underlying arable management or crop rotations include the ploughing-in of a perennial feed crop leading to large C inputs to the soil.

3.2.2 Nitrogen balance (NB)

In Fig. 3, the assessment of the distribution of the NB with the in- and outfluxes is presented. The averaged nitrogen influx (represented by the uncertainty ensemble mean) per hectare was estimated to be 212.3 ± 9.1 kg N ha⁻¹ yr⁻¹, with an IQR from 203.3 to 220.0 kg N ha⁻¹ yr⁻¹, while nitrogen outfluxes were estimated, on average, to be 198.3 ± 11.2 kg N ha⁻¹ yr⁻¹, with an IQR from 191.4 to 204.0 kg N ha⁻¹ yr⁻¹ (Fig. 3), leading to a net N accumulation in the soil of 14.0 ± 2.1 kg N ha⁻¹ yr⁻¹, with an IQR from 11.9 to 16.0 kg N ha⁻¹ yr¹.

Figure 3 Nitrogen balance for cropland cultivation for the region of Thessaly: (a) total NB (in k t N) and (b) averaged NB (in kg N ha⁻¹). Data are averaged for the years 2012-2016. Horizontal lines indicate the mean (red), median, and minimum and maximum of the distribution.

Table 5 Nitrogen balance (per hectare). Summary of the assessment and uncertainty analysis of the NB fluxes (per hectare) of cropland cultivation in the region of Thessaly, Greece. ¹ N inputs as the sum of the absolute values of all input fluxes of the 500 simulations. ² N outputs as the sum of the absolute values of all the output fluxes of the 500 simulations. ³ N stock changes as the difference between the input and output fluxes of each of the 500 simulations. ⁴ Gaseous emissions are the sum of N₂O, NO, N₂ and NH₃ fluxes. ⁵ Aquatic flux is nitrate leaching (NO${}_{\mathrm{3}}^{-}$).

The N influx was composed of the input of synthetic fertilizer of 80.2 ± 4.8 kg N ha⁻¹ yr⁻¹ (IQR 76.6 to 82.7) and the input of organic fertilizer of 80.2 ± 3.6 kg N ha⁻¹ yr⁻¹ (IQR from 77.5 to 82.7), followed by the biological nitrogen fixation (BNF) via legumes, estimated to be 45.6 ± 4.3 kg N ha⁻¹ yr⁻¹ (IQR from 43.7 to 47.7), and nitrogen deposition of 6.3 ± 0.8 kg N ha⁻¹ yr⁻¹ (IQR from 6.0 to 6.8). Thus, almost 75 % of the nitrogen input influx is related to the fertilization (mineral and organic), whilst the minor part that corresponds to nitrogen fixation and deposition approximates to 25 %.

The N outflux consist of gaseous N fluxes of 55.4 ± 8.8 kg N ha⁻¹ yr⁻¹ (IQR from 48.9 to 61.6); aquatic N fluxes (only nitrate leaching into surface waters was considered) of 14.1 ± 4.5 kg N ha⁻¹ yr⁻¹ (IQR from 11.0 to 17.0); and N fluxes of N removed in yields, straw and feed of 128.8 ± 8.5 kg N ha⁻¹ yr⁻¹ (IQR of 125.2 to 131.7) (see Fig. 4 and Table 5). Based on the aforementioned results, all gaseous and aquatic N fluxes correspond to about 28 % and 7 % of the N output flux respectively, while the far-largest N output flux was N removed in yields, straw and feed, representing almost 65 % of the N outflux (Fig. 4).

Figure 4 Averaged annual nitrogen balance (inner ring of the pie diagram) and their decomposition into the various components of the N fluxes (outer ring of the pie diagram); all data are summarized in Table 5.

The simulated gaseous fluxes were composed of N₂O emissions estimated to be 2.6 ± 0.8 kg N₂O–N ha⁻¹ yr⁻¹ (IQR from 2.1 to 3.1), NO emissions of 3.2 ± 1.5 kg NO–N ha⁻¹ yr⁻¹ (IQR from 2.0 to 4.1), N₂ emissions of 15.5 ± 7.0 kg N₂–N ha⁻¹ yr⁻¹ (IQR range from 9.9 to 20.7) and NH₃ emissions of 34.0 ± 6.7 kg NH₃–N ha⁻¹ yr⁻¹ (IQR from 29.3 to 36.9). Ammonia volatilization represents the largest share (61.48 %) of gaseous N losses, with the highest densities in the emission distribution between approximately 25 and 35 kg N ha⁻¹, followed by di-nitrogen losses (28.03 %) of gaseous N losses, with a much wider emission variability in the distribution, followed by NO₃ (5.79 %) and N₂O (4.7 %). Figure 5 shows the overall NB in a waterfall diagram adding up cumulatively all in- and outfluxes, illustrating the uncertainty distribution of each flux contribution. The waterfall diagram illustrates the overall outcome of the NB, an N accumulation into the soil as the difference between all outfluxes minus all influxes.

Figure 5 Waterfall representation of the result distributions of the different nitrogen in- and outfluxes of the cropland cultivation in Thessaly. Vertical lines in the distributions indicate mean values of the corresponding N flux. Red indicates gaseous outfluxes; blue indicates aquatic fluxes; green indicates biomass yield and feed production outfluxes; and brown indicates indicates N influxes such as synthetic N fertilizer, N manure, biological N fixation (BNF) and N deposition. The resulting N sink of the nitrogen balance (based on distribution means) is −13.8 kg N ha⁻¹ yr⁻¹. A negative value indicates flux into the soil.

The nitrate-leaching mean estimate was 14.1 ± 4.5 kg NO₃–N ha⁻¹ yr⁻¹ (IQR from 11.0 to 17.0) with a bell-shaped distribution.

Total yield and biomass (straw and feed) N export fluxes were 62.4 ± 4.4 kg N ha⁻¹ yr⁻¹, with uncertainty ranges from 59.9 to 65.1, consisting of yield N exports (grains and cotton balls) of 30.3 ± 1.7 kg N ha⁻¹ yr⁻¹ (IQR from 29.6 to 30.9) and straw and feed N exports of 36.1 ± 6.0 kg N ha⁻¹ yr⁻¹ (IQR from 34.9 to 37.6). The result distributions for yield N are notably bell shaped; for feed biomass N, they are very moderately bell shaped and well distributed within the bounds; and for straw N, they are very sharp within a comparably small interval.

Figure 5 illustrates the cumulative nitrogen fluxes of which the NB is composed as a waterfall diagram considering the mean of each component. The NB results in a net N sink of 13.8 kg N ha⁻¹ yr⁻¹ (see result distribution in Fig. 6) for the region, corresponding to an annual carbon sequestration of approximately 0.5 t C ha⁻¹ yr⁻¹ (as depicted in Fig. 2b) (see also the annual dynamics of the topsoil (3 cm) soil organic carbon and nitrogen distributions in Fig. 8).

Figure 6 Distribution of the overall nitrogen balance of the cropland cultivation in Thessaly: statistical analysis across all 500 individual NB results of the inventory simulations (mean 13.8 kg N ha⁻¹ yr⁻¹, median 13.7 kg N ha⁻¹ yr⁻¹) corresponding to the carbon balance in Fig. 2.

Figures 7 and 8 show the dynamics of the annual distribution of the gaseous and aquatic outfluxes, as well as the dynamics of the annual distributions of the topsoil (30 cm) soil organic carbon and nitrogen pools for the evaluation period 2011–2016.

Figure 7 Annual dynamics of the uncertainty distributions of the gaseous (a–d) and aquatic (e) N outfluxes for the period 2011–2016. Uncertainty bandwidth (blue band) is defined as the range between the Q25 and the Q75 quartiles; the green band (Q10 to Q90 interval) indicates the variance of the fluxes, neglecting the outliers of the distribution.

Figure 8 Annual dynamics of the uncertainty distributions of the soil carbon (a) and soil organic nitrogen (b) and the corresponding dynamics of the uncertainty distributions of the annual change rates of the total soil carbon and nitrogen pools (c, d) respectively.

4 Discussion

In this study, following the recommendation of Grosz et al. (2023), an assessment of the combined full C and N balances of a regional cropland agroecosystem is reported for the first time using inventory simulations with a process-based ecosystem model. The additional quantification of the associated modelling uncertainty of the balance simulations increases the trustworthiness of the study.

Up to present, process-based modelling studies have mainly focused on single-site applications, e.g. Daycent (del Grosso et al., 2005; Gurung et al., 2020), APSIM (Vogeler et al., 2013), CERES-EGC (Gabrielle et al., 2006; Dambreville et al., 2008; Heinen, 2006; Hénault et al., 2005), CERES-Wheat (Mavromatis, 2016), DNDC (Li, 2000) and LandscapeDNDC (Haas et al., 2013; Zhang et al., 2015; Klatt et al., 2015; Molina-Herrera et al., 2016). Fewer studies deploy models on the regional to national scale (del Grosso et al., 2005; Klatt et al., 2015; Kim et al., 2015) or the continental to global scale (del Grosso et al., 2009; Thompson et al., 2019; Franke et al., 2020; Smerald et al., 2022; Jägermeyr et al., 2021). All of these studies are subject to criticism, as stated by Grosz et al. (2023), as they are reporting, in general, on only one specific component or only a few components of the carbon or nitrogen cycles, such as soil carbon stocks or N₂O emissions, lacking any information on the full C and N balance.

There are very few cases where an attempt at a regional estimation of the NB has been made. The study by Schroeck et al. (2019) is the only previous attempt fulfilling the requirements of Grosz et al. (2023) in reporting on the full NB for a large alpine watershed in the Austrian Alps, characterized by arable production in the low-lying areas and grassland in the mountains, using a process-based model. In addition, Lee et al. (2020) tried to estimate nitrogen balances in Switzerland by alternating the cropping systems or management practices. There were also cases where the regional NB was estimated with the use of nitrogen balance equations (He et al., 2018). Recently, Zistl-Schlingmann et al. (2020) assessed the full N balance of alpine grasslands using the ¹⁵N stable isotope techniques.

In order to achieve a more concrete and complete analysis of the CB and NB that could be used for future policy development, an uncertainty analysis is considered to be necessary and mandatory. The IPCC guidelines demand UNFCCC (United Nations Framework Convention on Climate Change) reporting of the uncertainty quantification of any reported inventory study (IPCC, 2019). Recent publications have reported the deployment of different methods to assess and quantify the various sources of uncertainty in ecosystem modelling. Klatt et al. (2015) published a study on the impact of parameter uncertainty on N₂O emissions and NO₃ leaching at a regional scale. Houska et al. (2017) deployed the GLUE method (generalized likelihood uncertainty estimation) for the LandscapeDNDC model on a grassland site; other studies such as Lehuger et al. (2009a), Li et al. (2015) and Myrgiotis et al. (2018a) used the Bayesian model calibration and uncertainty assessment approach, which has been used in the current study as well.

4.1 Yield and feed production

LandscapeDNDC was validated in a study by Molina-Herrera et al. (2016) on cropland and grassland sites across Europe, showing good agreement in terms of reproducing observed above ground biomass and yield estimates. Similar model performance for the cultivation of commodity crops was reported by Kasper et al. (2019), Klatt et al. (2015), Molina-Herrera et al. (2017) and R. J. Petersen et al. (2021).

Lyra and Loukas (2021) used the REPIC model to estimate the crop growth and/or yield production of several crops in the Basin of Almyros, Thessaly. The simulated results were approximately 11 t ha⁻¹ clover, 3.3/3.5 t ha⁻¹ cereals and/or wheat, 3.8 t ha⁻¹ cotton, and 9 t ha⁻¹ maize; these values compare well to the results of our research, as shown in Table 3. The simulated results presented in our study are in line with the results by Voloudakis et al. (2015), simulating cotton production in seven different areas of Greece through application of the AquaCrop model. Similar results were reported by Tsakmakis et al. (2019).

There are few cases in the literature concerning yield simulations on a European level. Based on the yield datasets of FAO and EUROSTAT, Ciais et al. (2010a) estimated mean crop yields for the period 1990–1999 at the scale of EU-25 as 6.1 (FAO) and 5.3 (EUROSTAT) t DM ha⁻¹ yr⁻¹ respectively, which corresponds well to the results of our study. Haas et al. (2022) estimated with a model ensemble mean for crop yields for EU-27 a value of 4.41 ± 1.85 t DM ha⁻¹ yr⁻¹ for the period 1990–1999. Lugato et al. (2018) estimated cropland yield projections of 4.34 t DM ha⁻¹ yr⁻¹ (mean), ranging from 3.69 to 4.90 t DM ha⁻¹ yr⁻¹, with the DayCent model for EU-27, comparable to the 6.18 t DM ha⁻¹ yr⁻¹ average simulated crop yields of this study. The simulated yields in the current study vary from 3.3 to 9.9 t DM ha⁻¹ yr⁻¹ for the cases of cotton and maize respectively.

Higher yield estimates for the region of Thessaly in this study are certainly due to the inclusion of the legume feed crops in the rotations. This argument is supported by a recent meta-analysis by Lu (2020) that concluded that, on average, yield increases of 5.0 % to 25 % can be expected for various conditions if residues are completely returned to the field as compared to no-residue-return systems. Similar results were reported by Fuchs et al. (2020) and Barneze et al. (2020).

Following the recommendations of Grosz et al. (2023), our study has reported transparently all the major C and N fluxes of the region as simulated by the model. In our study, we have not calibrated the model against any observations; therefore, all simulation results will be discussed in comparison to other modelling results available. Presently, there is only one comparable modelling study available in the literature reporting and discussing the total N balance of a site or region, which we have used to compare our N balance to.

4.2 Carbon and nitrogen balance

Full N balance

At present, the studies of Schroeck et al. (2019) and Lee et al. (2020) are the only ones to be found on the Web of Science when searching the keywords “nitrogen”, “balance”, “process”, “based” and “modelling” simultaneously; these studies report a compilation of the nitrogen balance and all associated N fluxes for a site or region through applying a process-based ecosystem model, as demanded by Gosz et al. (2023).

Leip et al. (2011) reported the first nitrogen balance for Europe following a mixed approach combining the CAPRI (Common Agricultural Policy Regionalised Impact) model (a global economic model for agriculture) with different approaches estimating various nitrogen fluxes in arable land cultivation, but the approach lacks the explicit quantification of the different gaseous N fluxes. The study of Schroeck et al. (2019) overcame this hurdle and applied the process-based ecosystem model LandscapeDNDC to estimate the full regional nitrogen budgets, including all the fluxes of different ecosystems (cropland, grassland and pastures) and climatic zones, of a watershed in Austria. That has been the first attempt at estimating and reporting all the N fluxes possible, as demanded by Gosz et al. (2023).

The N balance estimate in Schroeck et al. (2019) for a catchment in Austria and the N balance reported in our study compare very well despite the inherent differences in land management and N inputs. As highlighted by Grosz et al. (2023), such intercomparisons demonstrate the different model behaviours when applied to different ecosystem. In our study, we see the partitioning of the N outfluxes from our arable system in similar shares as reported by Schroeck et al. (2019) for the arable land.

The N₂O estimates in Schroeck et al. (2019) and in the current study are of a comparable level. We estimated N₂O emissions of 2.6 kg N ha⁻¹ yr⁻¹, while Schroeck et al. (2019) report 1.51 kg N ha⁻¹ yr⁻¹, which is about 40 % lower. The NO fluxes differ significantly since we reported a mean value of 3.2 kg NO–N ha⁻¹ yr⁻¹, while Schroeck et al. (2019) reported 0.08 kg NO-N ha⁻¹ yr⁻¹. On the one hand, this is related to some recent model advances which were made during this study and which elevated the NO production in LandscapeDNDC (Molina-Herrera et al., 2017), and on the other hand, this is due to the high share of organic N fertilization in our study fostering NO emissions. Ammonia volatilization differs substantially between the two studies – while our study reports 34 kg NH₃–N ha⁻¹ yr⁻¹, Schroeck et al. (2019) reported moderate emissions of 0.23 kg NH₃–N ha⁻¹ yr⁻¹. The strong NH₃ volatilization in our study is mostly driven by the high pH values of the soils in the region of Thessaly (pH values from 6.5 to 8.2, with a considerable spatial variation; Greek Soil Map: https://iris.gov.gr/SoilServices/js/pdf/SOIL%20MAP%20OF%20GREECE%20e-SOILBOOK.pdf, last access: 17 April 2019, 2015) and the comparably high manure inputs into the arable system in our study, while in the research of Schroeck et al. (2019), the manure was preferably applied only to the grassland systems, and mineral fertilizers were applied to the arable land. Concerning the NO₃, Schroeck et al. (2019) reported a value of 45.3 kg NO₃–N ha⁻¹ yr⁻¹, which is 3 times higher than that estimated in this study (14.1 kg N ha⁻¹ yr⁻¹) considering N inputs of approximately 160 and 212.3 kg N ha⁻¹ yr⁻¹ respectively. Even though 50 % of the arable land in our study was irrigated, the resulting water percolation rates in our study were, by far, lower than the percolation simulated in the study of Schroeck et al. (2019) as the Austrian pre-alpine catchment received nearly double the annual precipitation.

The N balance modelling study of Lee et al. (2020) made estimations of a national cropland N balance for Switzerland using an upscaling method based on process-based site simulations with the DayCent model, distinguishing the management information of the considered cropping systems, e.g. fertilizer rates, tillage or land cover change. The study reported for conventional cultivations (averaged across 20 years) the yield-related N outfluxes, NO₃ leaching and gaseous N emissions, accounting for about 60 %, 36.1 % and 4.1 % respectively of the total N outputs. Lee et al. (2020) did not report the different gaseous N fluxes even though the DayCent model must have simulated all of them. Although the yield-related N outflux is in accordance with our result of 64.95 %, there seems to be a discrepancy in the reported gaseous and aquatic N flux contributions as we report 27.94 % for gaseous and 7.11 % for NO₃ leaching in our study. As demanded by Gosz et al (2023), we can elaborate upon different preferences in simulated N outflux partitioning (36 % NO₃ and 4 % gaseous losses for DayCent versus 7 % NO₃ and 28 % gaseous losses for LandscapeDNDC) due to the different simulation models and regionalization and upscaling approaches, as well as the different soil, climatic and management conditions included in the respective studies.

Velthof et al. (2009) used the MITTERA-EUROPE model and method, based on the concoction of the GAINS and CAPRI models, to estimate N fluxes of European soils on a NUTS2 scale, making use of European datasets and literature coefficients, where the fertilizer application and management were similar to our methodology. The average N input–output balance was calculated as 117 kg N ha⁻¹ yr⁻¹, composed of manure at 49 kg N ha⁻¹ yr⁻¹, synthetic fertilizer at 58 kg N ha⁻¹ yr⁻¹ (in the current study, for both cases equating to 80.2 kg N ha⁻¹ yr⁻¹), biological nitrogen fixation at 2 kg N ha⁻¹ yr⁻¹ (our research indicates 45.6 kg N ha⁻¹ yr⁻¹) and N deposition at 7 kg N ha⁻¹ (current study indicates 6.3 kg N ha⁻¹ yr⁻¹). In contrast to our study, the reported output fluxes are 8 kg NH₃–N ha⁻¹ yr⁻¹ for NH₃, 2 kg N₂O–N ha⁻¹ yr⁻¹ for N₂O, 2 kg NO_x–N ha⁻¹ yr⁻¹ for NO_x, 51 kg N₂–N ha⁻¹ yr⁻¹ for N₂ and 7 kg NO₃–N ha⁻¹ yr⁻¹ for NO₃ leaching, while the results presented in our study are 34.0 kg NH₃–N ha⁻¹ yr⁻¹ for NH₃, 2.6 kg N₂O–N ha⁻¹ yr⁻¹ for N₂O, 3.2 kg NO_x–N ha⁻¹ yr⁻¹ for NO_x, 15.5 kg N₂–N ha⁻¹ yr⁻¹ for N₂ and 14.1 kg NO₃–N ha⁻¹ yr⁻¹ for NO₃ leaching. Additionally, the yield output is estimated as 48 kg N ha⁻¹ yr⁻¹. Again, we see a different preference in N outflux partitioning for large shares in gaseous N fluxes versus small NO₃-leaching shares; the difference compared to the results presented in our study is related to the different input data used for initialization and driving of the model based on regional statistics and the use of a biogeochemical model versus emission factor approaches.

He et al. (2018) assessed the soil N balance for a time span between 1984 to 2014 based on the N budget equations (N input – N output) using multiple coefficients from the literature in order to estimate the nitrogen input and output fluxes of six grouped regions in China. The used datasets were acquired from national authorities and include cropping land and yields, synthetic fertilizers, animal heads, soil types, etc. The N synthetic fertilizer input is, on average, 182.4 kg N ha⁻¹; the organic fertilizer is 97.3 kg N ha⁻¹; N fixation is estimated as 16.8 kg N ha⁻¹; and the atmospheric deposition is estimated as 22 kg N ha⁻¹. Almost half of the total averaged N output losses, 48.9 %, was attributed to crop uptake, while the respective gaseous losses were 19.9 % for N₂, 17.3 % for volatilized NH₃, 1.2 % for N₂O and 0.7 % for NO. As for the NO₃ leaching, the share was 5.8 % of the total output N fluxes. These reported N outflux proportions were notably comparable to our study. The differences in the N uptake data remain and are mainly due to the differences in the crops and management.

As reported by the OECD (OECD, 2020), the net averaged nitrogen balance of the area of our study is 11.6 kg N ha⁻¹ yr⁻¹ input to the soil, which corresponds very well to the simulated mean nitrogen balance as an influx of 13.8 kg N ha⁻¹ yr⁻¹ (IQR 11.9 to 16.0) into the soil.

So far, the discussion of the presented N balance and N outfluxes compares well to most of the available studies reporting N balances; however, one modelling study reported different N outflux partitioning between gaseous and NO₃-leaching fluxes. For a more detailed intercomparison of the overall quality of our C and N fluxes, we aim to compare our results to various studies, addressing individual components of the C and N balance and associated fluxes.

SOC stocks

Haas et al. (2022) reported results of a European inventory simulation of soil carbon stocks and N₂O emissions using a model ensemble. The study was deployed in a baseline simulation across EU-27 with a similar residue management as compared to our study, resulting in very stable carbon stock dynamics over a long period (1950–2100). In this study, the estimated carbon sequestration of 0.5 (UA mean and median) ± 0.3 t C ha⁻¹ yr⁻¹ is mainly caused by the inclusion of legume feed crops within the crop rotation, leading to increased litter production and C input into the soil (K. Petersen et al., 2021; Fuchs et al., 2020; Barneze et al., 2020). Haas et al. (2022) reported a management scenario with 100 % of the crop litter remaining on the field, leading to averaged C sequestration rates of over 1 t C ha⁻¹ yr⁻¹ across EU-27. As the residue management in this study is between the baseline and buried scenario of Haas et al. (2022), our results compare well to the results reported in that study.

Other modelling studies such as Lugato et al. (2014) reported C sequestration rates for the conversion of cropland into grassland ranging between 0.4 and 0.8 t C ha⁻¹ yr⁻¹. Lugato et al. (2014) reported averaged SOC change rates for a cereal straw incorporation scenario for EU-27 of 0.1 t C ha⁻¹ yr⁻¹ (estimates from 2000–2020).

The SOC dynamics reported in this study show a stable carbon dynamic in the soil within the simulation time span (2009–2014) with only 3 years of model spin-up. The initialization of the various carbon pools with the SOC data from the soil database is balanced by the average litter production of the deployed crop rotations. The SOC increase in 2015 and 2016 is due to climatic conditions and higher litter inputs simulated by the model.

N₂O emissions

This study reported estimates of N₂O emissions of 2.6 ± 0.8 kg N₂O–N ha⁻¹ yr⁻¹ (IQR from 2.1 to 3.1) for a mixed-crop and/or legume feed crop rotation, which were well above the estimates resulting from IPCC tier-I direct-emission factors; IPCC estimated 1.6 kg N₂O–N ha⁻¹ yr⁻¹ when applying approximately 160 kg N ha⁻¹ yr⁻¹. The higher N₂O emission strength of the modelling is likely to result from emission peaks after irrigation due to low anaerobicity (Grosz et al., 2023; Janz et al., 2022). Cayuela et al. (2017) conducted a meta-analysis of the direct N₂O emissions for a number of cropping systems for a Mediterranean climate, where the emission factors (EFs) were altered under different fertilization and irrigation conditions. Higher fertilization rates led to higher EFs (0.82 % less than the 1 % of IPCC). Additionally, irrigated and intensively cultivated crops had higher EFs than rain-fed crops (up to 0.91 % depending on the irrigation method). The relatively high EF of maize in this study could possibly be attributed to the irrigation without the application of water-saving methods and the, on average, higher N application rates .

The LandscapeDNDC validation study of Molina-Herrera et al. (2016) reported for the Italian site Borgo Cioffi (Mediterranean climate; Ranucci et al. (2011)) annual N₂O emissions of 2.49 kg N₂O–N ha⁻¹ yr⁻¹, while two sites in southern France showed annual N₂O emissions from 0.52 to 3.34 kg N₂O–N ha⁻¹ yr⁻¹. The N₂O emission estimates of our study were higher than the results reported by Haas et al. (2022) using a multi-model ensemble estimating average soil N₂O emissions from European (EU-27) cropping systems for the period 1980–1999 (1.46 ± 1.30 kg N₂O–N ha⁻¹ yr⁻¹) under conventional (baseline) management and comparably average N input. Klatt et al. (2015) reported for a regional inventory (Saxony, Germany) a mean N₂O emission of 1.43 ± 1.25 kg N₂O–N ha⁻¹ yr⁻¹.

Overall, the reported N₂O flux component of our study compares well to the findings reported in the literature. As criticized by Grosz et al. (2023), many studies only focus on the performance of the models in simulating N₂O emissions, and these models were even calibrated for this purpose. Without reporting all the other N fluxes from the models, this focus on and calibration for only one quantity can easily lead to inaccuracies with regard to other components of the N cycle as they may not be checked for consistency (Janz et al., 2022).

Nitrate leaching

This study reported average NO₃-leaching fluxes (only nitrate leaching into surface waters) of 14.1 ± 4.5 kg NO₃–N ha⁻¹ yr⁻¹. Reported nitrate-leaching observations for the region or Greece generally could not be found in the literature; de Vries et al. (2011) estimated the NO₃ leaching with the use of four different models, with varying values from 5 to 40 kg NO₃–N ha⁻¹ yr⁻¹ for the area of our study. These high values could be explained by the fact that they correspond to both groundwater and runoff. Molina-Herrera et al. (2016) reported for the LandscapeDNDC validation study cropland nitrate-leaching fluxes of approximately 7 to 88 kg NO₃–N ha⁻¹ yr⁻¹. In addition, in the research of Molina-Herrera et al. (2017), the described NO₃-leaching results varied from 13 to 8 kg NO₃–N ha⁻¹ yr⁻¹, showing higher values with regard to the precipitation and fertigation. The most comparable site, Borgo Cioffi, resulted in a comparable annual NO₃-leaching flux of 18.62 kg NO₃–N ha⁻¹ yr⁻¹.

Klatt et al. (2015) reported in an uncertainty assessment for a regional inventory (Saxony, Germany) leaching rates of 29.32 ± 9.97 kg NO₃–N ha⁻¹ yr⁻¹ for a wheat–barley–rapeseed rotation simulated by the LandscapeDNDC model. The agricultural system and management regime are comparable; higher NO₃-leaching rates were most likely due to high N fertilization rates in combination with higher annual precipitation in the region, leading to more intense percolation and therefore to stronger leaching of available NO₃, while, in our study, the fertilization regime was more lean, such that soil nutrient competition was higher, and available nitrate was more likely to be immobilized by plant uptake. Myrgiotis et al. (2019) reported in a similar assessment an NO₃-leaching factor (LF) mean of 14 % (± 7  %) for their region; in comparison, we report a mean NO₃-leaching factor of 7 %.

NO emissions

In the current study, the model-estimated NO emissions were, on average, 3.2 ± 1.5 kg NO–N ha⁻¹ yr⁻¹. Butterbach-Bahl et al. (2009) performed the very first European inventory of soil NO emissions using a modified version of DNDC, reporting low NO emission rates mostly below 2 kg NO–N ha⁻¹ yr⁻¹. Molina-Herrera et al. (2017) recently reported a full NO emission inventory for Saxony, Germany, compiling annual NO emissions from agricultural soils ranging from 0.19 to 6.7 kg NO–N ha⁻¹ yr⁻¹, as simulated by LandscapeDNDC. The study reported the model performance in simulating soil NO emissions at more than 20 different sites. The study of Schroeck et al. (2019) reported for a regional inventory of arable soils in Austria, as simulated by LandscapeDNDC, annual NO emissions of 1.0–1.5 kg NO–N ha⁻¹ (for the year 2000), while empirical approaches such as that of Stehfest and Bouwman (2006) estimated emissions of similar magnitude. Zhang et al. (2015) reported in a model intercomparison and validation study (including three ecosystem models) of NO and N₂O fluxes consistent simulation results for the LandscapeDNDC model, with NO emission strengths of cropland soils between 1 and 3 kg NO–N ha⁻¹ yr⁻¹ across the sites.

NH₃ emissions

Schroeck et al. (2019) stated that validation studies of NH₃ volatilization for any biogeochemical model are very rarely reported in the literature, mainly due to the complexity and a lack of flux observations at high spatial and temporal resolutions.

In our study, we estimate soil NH₃ emissions of 34.0 ± 6.7 kg NH₃–N ha⁻¹ yr⁻¹. High NH₃ volatilization and emission rates can be explained by the predominating neutral to basal soil conditions (pH values of 7 and above) in the study region, favouring the Henry NH₄–NH₃ equilibrium towards higher NH₃ gases and enabling ammonia to diffuse out of the soil into the free atmosphere.

The IPCC emission factor (EF) method for NH₃ volatilization reports estimates of 20 % of N input into the soil being volatilized as NH₃. For our study, the IPCC methodology for NH₃ would lead to 32 kg NH₃–N ha⁻¹ yr⁻¹, which is well in line with the simulated result.

Sidiropoulos and Tsilingiridis (2009) estimated a national livestock-originated NH₃ emission corresponding to approximately 22 kg ha⁻¹ yr⁻¹ for the region of Thessaly.

There is a number of national NH₃ inventories which could be considered to be detailed and well-studied, such as the ones in Denmark, the Netherlands, Europe, the UK and the US. In Denmark, Geels et al. (2012) used the DAMOS model to estimate the Danish NH₃ emissions (crop, grass and manure manipulation); in the five regions under study, the values ranged from a very small quantity to 17.4 kg NH₃–N ha⁻¹ yr⁻¹.

As discussed by Sutton et al. (2013), the majority of the NH₃ emissions come about as a result of agricultural production and are considerably impacted by climate influence. In the case of NH₃ volatilization, it could almost double with every 5 °C in temperature given certain complex thermodynamic dissociation and solubility, whilst soil NH₃ emission is influenced by the available water quantity, allowing the NH_x dissolution and use by microbial organisms, which afterwards leads to decomposition.

4.3 Uncertainty analysis and quantification

Santabárbara (2019) used the MCMC algorithm to estimate the joint parameter distribution of the fundamental bio-geochemical process parameters in LandscapeDNDC when simulating soil C and N fluxes. Propagating these joint parameter distributions through the model (by sampling 500 joint parameter distributions and performing inventory simulations with each parameter set within the model) for estimating the regional C and N fluxes lead to various distributions for any model result on a regional scale. Statistical analysis calculating the mean, median and interquartile range (Q25 to Q75) determines the best estimates and the uncertainty range of any model output on a regional scale, demonstrating the superiority of the method for assessing any ecosystem response by modelling instead of reporting single results. This is a novel approach that, to our knowledge, has not been reported before in the literature for the full carbon and nitrogen balance and has not been applied to regional simulations by any process-based model.

In this study, the estimated UA mean and median of the carbon sequestration of 0.5 ± 0.3 t C ha⁻¹ yr⁻¹ are associated with an uncertainty range from 0.4 to 0.7 t C ha⁻¹ yr⁻¹, which compares well to the spatial uncertainty of C sequestration in the study of Haas et al. (2022). The approach used in this study enabled us to assess the carbon and nitrogen balance of the Bayesian calibration method used by Lehuger et al. (2009b) for the enhancement of the CERES-EGC model parameterization (reduction of the a priori parameter distribution), as well as the quantification of the uncertainty of the simulated N₂O emissions at different sites. The estimated fluxes of the different sites resulted in a range between 0.088 to 3.672 kg N₂O–N ha⁻¹ yr⁻¹, with values of 0.066 to 0.115 kg N₂O–N ha⁻¹ yr⁻¹ and 1.676 to 5.874 kg N₂O–N ha⁻¹ yr⁻¹ for the Q05 and Q95 quantiles respectively, with an averaged value of 1.04 kg N₂O–N ha⁻¹ yr⁻¹, which is lower than the result of the current study but still of the same order of magnitude.

Klatt et al. (2015) quantified a parameter-induced uncertainty analysis at a regional scale by applying the same process model for simulating N₂O emission and NO₃-leaching inventories as that in our study. The region was represented by 4000 polygons of arable land (Saxony, Germany) for crop rotations of barley, wheat and rapeseed, with differing climatic conditions. The results of Klatt et al. (2015) display a likelihood range of 50 % (the IQR range between Q25 and Q75) for N₂O emissions from 0.46 to 2.05 kg N₂O–N ha⁻¹ yr⁻¹, which is in good agreement with our results of 2.1 to 3.1 kg N₂O–N ha⁻¹ yr⁻¹. The average N₂O emissions are 1.43 kg N₂O–N ha⁻¹ yr⁻¹, which is comparable to the result of our study (mean: 2.6 kg N₂O–N ha⁻¹ yr⁻¹, median: 2.5 kg N₂O–N ha⁻¹ yr⁻¹ across approximately 1000 polygons). As for leached NO₃, Klatt et al. (2015) reported mean leaching rates of 29 kg NO₃–N ha⁻¹ yr⁻¹ (IQR from 24.5 to 36.0), which is higher compared to the results of our study (mean: 14.1 kg NO₃-N ha⁻¹ yr⁻¹, median: 13.6 kg NO₃-N ha⁻¹ yr⁻¹) (IQR from 11 to 17). Despite the difference in climatic and soil conditions, both uncertainty analysis studies reported similar regional estimates and uncertainty ranges for N₂O emissions and NO₃ leaching.

Butterbach-Bahl et al. (2022) reported the influence of management uncertainties on compiling national inventories of CH₄ and N₂O emissions for various rice cultivation systems in Vietnam. The study applied a sampling technique that varies model input data within a given range and analysed the influence on the assessed CH₄ and N₂O emission strengths. As the underlying cropland systems were fundamentally different, the assessed uncertainty ranges were comparable, and the study supports our approach of focusing on uncertainty ranges rather than single values.

5 Conclusions

In this research, we presented, for the first time, a regional inventory of the full carbon and nitrogen balance, including all sub-components of these fluxes simulated by a process-based model. Additionally, the study has fulfilled the demand of always reporting the associated uncertainties for any modelling results being published in literature. This supports the trustworthiness of the reported results for the C and N balances.

Comparing the modelled N balance with a similar approach that models the full N balance with all associated fluxes for a catchment in pre-alpine Austria leads to the conclusion that partitioning the N outflux into the different N flux components is more inherent to the LandscapeDNDC model itself, as used in both studies, than it is induced by the two very different agricultural and climatical systems. Nevertheless, specific N outfluxes between the two studies show large differences (e.g. NH₃ volatilization), which is purely caused by model processes due to different soil pH values. Comparing to a less granular and detailed study of the N balance for Switzerland gives a first impression of the differences to be expected in modelling the arable N balance with various different models. The discussion of such results will become more lively and maybe controversial as soon as more comparable studies using different models become available.

In addition, a full uncertainty analysis is presented based on the Metropolis–Hastings algorithm, where a parameter subset and input data perturbation were sampled and simulated, resulting in various probability density functions (PDFs) for each one of the N and C balance fluxes, building a full uncertainty analysis of the modelled results. This helps to build trustworthiness in modelling assessments and estimates of the balances and of the model behaviour.

As demanded by the nitrogen-modelling community, all of the above constitute the novelty of the conducted research that could be seen as a prototype for analysing and reporting N cycling in agro-ecosystems in the future.