Reviewer 1 (Emanuele Lugato)

Dear editor

The paper of Sifounakis et al., presents a regional simulation of C and N budget in Thessaly (Greece), with the widely used biogeochemical model landscapeDNDC. The strength of the paper is that is driven by detailed regional information and uses an uncertainty analysis based on the Markov Chain Monte Carlo (MCMC) Metropolis–Hastings algorithm. Information about model implementation and parameterization, as well as uncertainty calculations, are quite useful for the GHG inventory community (and beyond), considering the attempt to move toward tier 3 approaches under the EU LULUCF regulation.

Up to the discussion, the paper is substantially well developed and the overall N fluxes (and C to a lesser extent) look reliable according to other studies and my personal experience.

I have, anyway, a series of questions and concerns that, I hope, can lead to some improvements of the manuscript.

• Methods

The European Soil Database (ESDB v2.0, 2004), used as input soil dataset, is quite old indeed. New gridded value are available at ESDAC or alternative at ISRIC (soil grid), which authors may think to consider for the next studies (if not for this).

This is indeed an important criticism about the regional input data used in the study. As most regional inventory studies consume 75% of the study effort in collecting, aggregating and preparing of regional input data for model initialization and to drive the model along the simulation time span. As the arable land management (crop rotations, residue management and fertilization/manuring) poses the biggest uncertainty for the N balance, we have focused to use the available resources in this study to use most recent regional information of arable land management for the region, and use an existing model initialization dataset from the EU project NitroEurope based on the European Soil Database (ESDB v2.0, 2004). To our understanding, the improvements in using a newer version of the Soil Database or using ISRIC soil grid data, would improve the accuracy of the inventory calculation to a lower extend than improvements in the description of the real arable land cultivation would have on the carbon and nitrogen cycling of the region. Nevertheless we aim to follow this comment and to test the influence of different sources of soil data for the simulation of the arable land carbon and nitrogen cycling.

While the uncertainty is an important and valuable part of the work, some clarifications and details are needed, considering that Authors refer to a previous work ‘Santabarbara (2019)’ missing in the reference list. In the text is stated: “In the current analysis, 500 joint parameter sets were sampled from the posterior distributions in combination with input data perturbations as reported by Santabarbara (2019)”.

Oh sorry for the missing reference and the invalid citation. We have added the reference of Santabarbara (2019) to the manuscript. We use the Mendeley citation software and we will check all citations and references again.



Is not totally clear if MCMC Metropolis–Hastings is performed for this study or if posterior distributions to generate the 500 joint drawing (runs) are derived from a previous work. I also assume that posterior distributions for the initial soil conditions and drivers are generated for each of 1000 polygons simulated (in total I count 1000*10*500 simulations). Moreover, the 24 most sensitive process parameters to gaseous N fluxes are not mentioned and would be interesting to report them in the supplementary.

We have added a paragraph starting in line 238 explaining the methodology in more details: In a previous study by Santabarbara (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 supplementary material), the most sensitive soil initial data (soil profile data on bulk density, soil organic carbon content, pH value) and the most sensitive management information (fertilization and manure N rates, tilling depth) to aquatic and gaseous N fluxes from arable soils. This was digested in the MCMC simulation sampling a combination of 24 parameter values, 3 values of soil initial data and 3 management information. The sampling of the soil initial data as well as 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 Santabarbara (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 / input data perturbation sets from the posterior distributions as reported by Santabarbara (2019) and we deployed them in simulations (propagation through the model) for the regional inventory leading to 500 inventory simulations.  

I generally recommend to improve this section providing more details, considering eventually a flowchart to facilitate the reader.

In the preparation of the manuscript we aimed to focus more on the deployment of the uncertainty analysis rather than details on the method. We appreciate your comment and will add a section in the supplementary material describing the details of the i) sensitivity analysis, ii) the theory and application of the MCMC simulation to generate the posterior distributions of parameter / disturbance sets and iii) the propagation of the posterior distributions through the model leading to result distributions.

If I well understood, each crop rotation (R1-R5) is simulated in 100% of the agricultural land and then each crop, in a given year, is weighted based on statistics about corresponding cultivated area. Therefore, for example (table 1-2), the clover increases from 0.15 to 0.39 from 2012 to 2013, but the additional area is implicitly coming totally from a preceding rotation with winter wheat.  That may lead to some approximations that can be acceptable, but could be also mentioned.

Yes, your understanding of the methodology is correct. It is well known from literature, that the agricultural performance as well as the strength of the soil carbon and nitrogen cycling strongly depends on the interaction with crop rotations. Crop rotations will return various nutrients in different quality to the soil and they are key to prevent / interrupt pest and disease cycles. While the former is not implemented in LandscapeDNDC, effects of crop rotations to improve soil health by filling the soil organic matter pools from root litter and above ground residues from different crops (depending on the different litter carbon to nitrogen ratios to be transferred the various soil carbon pools) are well represented in the model (see Haas et al., 2022 (STOTEN) Long term impact of residue management on soil organic carbon stocks and nitrous oxide emissions from European croplands, https://doi.org/10.1016/j.scitotenv.2022.154932). The optimal solution would be to use spatial high resolution detailed information (such as EU invekos data) on field scale to take rotation effects into account. Data protection rules and data availability constraints prevents this such that modelling efforts need to simplify. In our opinion, it we try with the construction of the five rotations to come close as possible to reality while keeping complexity low as possible. In contrary to our very complex approach, recent global and continental inventory simulations such as Jägermeyr, J., et al, 2021 (Climate impacts on global agriculture emerge earlier in new generation of climate and crop models. Nature Food, 2, no. 11, 873-885, doi:10.1038/s43016-021-00400-y. ) perform several single crop monoculture simulations over long time spans. They account only for crop residues from the same crop from previous years repetitively. These approaches are highly vulnerable to artificial soil carbon and nitrogen depletions and accumulations.  

In the review, we would add a paragraph explaining in more details the construction of the crop rotations and their interactions.

• Carbon budget

In general, I found some semantic problems in C budget nomenclature. GPP+ manure cannot be called ‘C input' to soil as half of GPP is respired autotrophically and a consistent part removed from field. Moreover, what is called carbon fluxes from the soil to the atmosphere (TER) is properly the ‘ecosystems respiration’, since it includes both autotrophic and heterotrophic respiration. The terminology should be carefully revised throughout the text.

The C budget shows an average SOC sequestration rate of 0.5 t/ha per year, which seems a quite high rate of accumulation to my experience. I’m wondering how this number is affected by the model initialization in term of i) uncertainty of initial values, especially if coming from ESDB and ii) SOC pools partition, which is not mentioned/described. Not having a long-term spin-up, the model could be far from its equilibrium state, leading to spurious trends.

The selected soil bio-geochemistry module in the LandscapeDNDC model works in contrast to e.g. the Daycent model different: i) The initial soil carbon value per soil layer in the input file is divided during initialization into the internal 4 carbon pools (depending on their decomposability). This process was well calibrated in the past. Ii) To ensure stability with respect to soc input data, the model adjusts the decomposition rates of the 3 soil carbon litter pools (within valid ranges) within the first 3 simulation years. This method has proven to ensure SOC stability and an equilibrium when simulating forest, grassland and arable systems. The model has been extensively tested agains the long term data from the Rothhamsted soil carbon experiments and a dataset from the Askov site in Denmark and a manuscript about details on soc stability and model equilibrium is under preparation.

The explanation that this SOC accrual is ‘mainly caused by the inclusion of legume feed crops within the crop rotation leading to increased litter production’ is not convincing, unless you assume that those fodder crops were not in the rotation before 2009.

The C sequestration rate of approx.. 0.5 t/ha per year reported in this study is higher than averages across literature for pure cropping systems. In our study the effect of the perennial legume feed crop within the rotation (which is on the field for 18 months) builds up vegetative carbon stored in roots of several tons per hectare. This carbon will be transferred into the SOC litter pools during the final harvest and tillage. The C sequestration without this grass effect diminish and overall be very depending on the residues management assumed. We have used residues return numbers from a recent study by Haas. et al. (2022) analyzing the effects of SOC dynamics under different residue management scenarios across EU-27.

Here we are faced with the overall uncertainty issue as the data sources differ. Yes it is right, that our soil initialization will most likely not recognize the legume crop shares from the local cropping statistics. Adjusting soil initial data was not an option as this will for sure introduce uncertainty.

• Figures

The quality and details of the figures should be checked with care. For instance, Fig.3 is repeated 3 time without sub-figure labeling (a,b) and, in Fig.8 , the legend is covering the lines.  The Fig. 4 is not so appealing and, maybe, redundant having Fig. 5.

We will carefully revise the illustrations, captions and sub-figure labeling.

Since the authors made a regional model application, I would have been interested to see some maps, rather than eventually temporal trends that could be accommodated in supplementary material. 

We are well aware of this point. As the study is already quite extensive (regional inventory simulation, presentation of the full C and N balance and presentation of the model uncertainty) we had to limit the analysis for this paper. For the authors, it has been a priority to report the full C and N balance with all fluxes. As there was only one modelling paper found reporting the full N balance (of a region) we are aware that many questions about the robustness of the modelling effort will arise and therefore we decided to report the uncertainty analysis. A detailed reporting and analysis of spatial results (spatial patterns of results and geospatial analysis of soil and climate drivers) would cover a substantial part of the paper, which we will address in a follow up paper for sure.

• Discussion

In my opinion, this is the part of the paper that can be improved more. Most of time, the discussion is developing around the comparison of simulated DNDC fluxes with other modelling or empirical-driven approaches done at different scale, time, in other continents or sites with different pedo-climatic conditions. We all know how fluxes, especially N, are dependent on the local context, therefore, I would only keep the comparisons with data on the specific region.  Eventually, all the other data reported by different studies for the different fluxes could be summarized in a graphical way, shortening the discussion and indicating the limits mentioned above.

We are aware that reporting N balances for a specific region is a challenge especially when no comparable studies neither experimental nor numerical studies are available. For Greece or the Mediterrian region, we did not find any comparable results. Neither did we find comparable literature reporting at least parts of the N cycle. Our objective was to present the N balance as a new standard when simulating and analyzing any part of the N cycle (Nitrate leaching, No2, No, NH3 emissiones, etc.). Therefore we have tried to present the analysis as robust and as precise as possible. This includes i) the presented N balance including all sub N fluxes compares well with the only other simulation study found (using LandscapeDNDC as well),  ii) the sub fluxes compare well with observations (even these sub fluxes belong to many different systems), iii) the inclusion of an extensive uncertainty analysis helps to ensure trustworthiness in the analysis.

I would have expected to understand more about the model behavior in space, which (and why) modifications have been introduced in this run and more about model sensitivity to the different parameters. For instance, a policymaker might be interested tp understand which input information (soil, management, crop etc.) should be improved more to reduce the uncertainty or, if uncertainty has the same magnitude regionally.         

This is a very important issue especially as more and more stakeholders as for easy answers for mitigation. We will aim to address this in a follow up study where we will focus on geospatial result analysis and try to include the mitigation option as well.

I think the paper has good basis, with a consistent model parameterization based on regional data and reliable simulated fluxes, especially for the N cycle. General improvements are needed, especially in the uncertainty methodology description and the discussion, not very informative and comparing often pears and apples. 

Finally we want to state that in Germany there has a six years research initiative about Denitrification just recently ended with an opinion paper (submitted) concluding that the scientific community still lacks insights into the overall N budgets in agricultural systems. Only very few recent papers started reporting and analyzing the overall N balance especially using stable isotope measurements to address denitrification and N2 losses.

But for the modelling community to our knowledge, any model targeting N2O would always report the overall N budget. The Denitrification initiative concluded that there were no N budgets reported in the past such that conclusions on the performance of the different models could not be derived. (Except the on study by Schroeck using the LandscapeDNDC model).

The conclusion of the paper was to motivate the modelling community to report the overall N balance and not only sub fluxes.

The spin-up to achieve equilibrium is very short and the SOC results in figure 8 let assume that the system might be not in equilibrium (the changes are not explained in the discussion). I am also wondering about the assumptions made for the residue management, as this will be a crucial impact on the SOC changes. While this aspect was discussed for other studies, there is no analysis for the here presented results.

The selected soil bio-geochemistry module in the LandscapeDNDC model works in contrast to e.g. the Daycent model different: i) The initial soil carbon value per soil layer in the input file is divided during initialization into the internal 4 carbon pools (depending on their decomposability). This process was well calibrated in the past. Ii) To ensure stability with respect to soc input data, the model adjusts the decomposition rates of the 3 soil carbon litter pools (within valid ranges) within the first 3 simulation years. This method has proven to ensure SOC stability and an equilibrium when simulating forest, grassland and arable systems. The model has been extensively tested against the long term data from the Rothamsted soil carbon experiments and a dataset from the Askov site in Denmark and a manuscript about details on soc stability and model equilibrium is under preparation.

The C sequestration rate of approx.. 0.5 t/ha per year reported in this study is higher than averages across literature for pure cropping systems. In our study the effect of the perennial legume feed crop within the rotation (which is on the field for 18 months) builds up vegetative carbon stored in roots of several tons per hectare. This carbon will be transferred into the SOC litter pools during the final harvest and tillage. The C sequestration without this grass effect diminish and overall be very depending on the residues management assumed. We have used residues return numbers from a recent study by Haas. et al. (2022) analyzing the effects of SOC dynamics under different residue management scenarios across EU-27. Residues return values were selected as in the “Basecase” scenario for Haas at all (2022) showing a stable SOC dynamics for a 3 model ensemble.

I am not sure, why climate change mitigation is listed in the keywords. This is not part of this study. Actually, this study missed the opportunity to analyse the different fluxes against each other, as there are not many studies analysing bot: C- and N-fluxes. This provides the opportunity to analyse, if SOC gains will be compensated by N-emissions. Is the spatial distribution relevant for this compensation? ….

The use of mitigation is unfortune as the study is nor focusing on mitigation. The keyword was used to indicate the results to be used as a baseline for the N balance. We will correct this.

We are well aware of this point. As the study is already quite extensive (regional inventory simulation, presentation of the full C and N balance and presentation of the model uncertainty) we had to limit the analysis for this paper. For the authors, it has been a priority to report the full C and N balance with all fluxes. As there was only one modelling paper found reporting the full N balance (of a region) we are aware that many questions about the robustness of the modelling effort will arise and therefore we decided to report the uncertainty analysis.

A detailed reporting and analysis of spatial results (spatial patterns of results and geospatial analysis of soil and climate drivers) would cover a substantial part of the paper, which we will address in a follow up paper for sure.

I suggest to analyse the own results more detailed and provide an appropriate discussion of the results.

Finally we want emphasize the overall motivation to design the study and the manuscript in the way it has been presented submitted.

In Germany there has a six years research initiative founded by the German Science Foundation about Denitrification just recently ended with an opinion paper (submitted) concluding that the scientific community still lacks insights into the overall N budgets in agricultural systems. Only very few recent papers started reporting and analyzing the overall N balance (including all N sub fluxes) especially using stable isotope measurements to address denitrification and N2 losses.

But for the modelling community to our knowledge, any model targeting N2O would always report the overall N budget. No study so far did as the scientific community i) either asked for validation which could not be delivered due to lack of observations available, or ii) the model performance for the overall N balance or other N sub fluxes were not always as “good” as expected and authors feared criticism for their scientific work and findings. The Denitrification initiative concluded that there were no N budgets reported in the past such that conclusions on the overall performance of the different bio-geochemical models could not be derived! (Except the on study by Schroeck using the LandscapeDNDC model).

The conclusion of the paper was to motivate the modelling community to report the overall N balance and not only sub fluxes even some sub results may be assigned with high uncertainties. The community needs these results to identify shortcomings and derive decisions where to focus model adaptations / developments.

Some additional comments:

Abstract

There are no comments or results on the uncertainty analysis in the abstract. Lines 24-34 are only a list of results, but do not summarise the paper or aggregate the outcome of the study.

 We will revise this in the review.

Graphical abstract

The picture is nice, but it shows only the nitrogen fluxes, but not any carbon fluxes or any information about the uncertainty analysis. This means, that this does not summarise the study.

For us the novelty and therefore the focus was on the N budget. The C budget has been extensively reported by other studies before.

Methods

Table 1 indicates corn as maize, summarising food corn and silage maize. For modelling, however, there is a crucial difference, as the residue treatment differs and will affect large differences for the carbon input to the soil.

Unfortunately, the statistics for the region do not distinguishes between food corn and silage maize. We are aware that especially for the resiude management this is crucial. We did mention this in the paper and will  extend it in the review.

Table 2 suggest that the rotation changed, which does not make sense. Does this represent the crop coverage for the corresponding year or is a new rotation introduced each year? I assume this is due to the fact that there is a five year rotation applied on a 8 year simulation period with different representations. Is this possibly affecting the results?

It is well known from literature, that the agricultural performance as well as the strength of the soil carbon and nitrogen cycling strongly depends on the interaction with crop rotations. Crop rotations will return various nutrients in different quality to the soil and they are key to prevent / interrupt pest and disease cycles. While the former is not implemented in LandscapeDNDC, effects of crop rotations to improve soil health by filling the soil organic matter pools from root litter and above ground residues from different crops (depending on the different litter carbon to nitrogen ratios to be transferred the various soil carbon pools) are well represented in the model (see Haas et al., 2022 (STOTEN). Long term impact of residue management on soil organic carbon stocks and nitrous oxide emissions from European croplands, https://doi.org/10.1016/j.scitotenv.2022.154932).

We did not have any detailed information on crop rotations for the region. The only data we could derive were crop cultivation statistics. Therefore we had to derive a suitable crop rotation. At the same time, we had data of the share of arable land used per crop in the region. Therefore we duplicated the crop rotation in such a way (shifted them by one year) that each crop occurs only once per year in one of the five rotations. Then the final five rotations were extrapolated into the future.  

The optimal solution would be to use spatial high resolution detailed information (such as EU invekos data) on field scale to take rotation effects into account. Data protection rules and data availability constraints prevents this such that modelling efforts need to simplify. In our opinion, it we try with the construction of the five rotations to come close as possible to reality while keeping complexity low as possible. In contrary to our very complex approach, recent global and continental inventory simulations such as Jägermeyr, J., et al, 2021 (Climate impacts on global agriculture emerge earlier in new generation of climate and crop models. Nature Food, 2, no. 11, 873-885, doi:10.1038/s43016-021-00400-y. ) perform several single crop monoculture simulations over long time spans. They account only for crop residues from the same crop from previous years repetitively. These approaches are highly vulnerable to artificial soil carbon and nitrogen depletions and accumulations.  

Yes crop rotations will influence overall results, but we have checked with local farmer advisers for the most suitable rotation out of our 5 crops. Therefore will test the sensitivity of the influence of the construction of the crop rotation in the review for the paper and e.g. add this to the supplementary material.

In the review, we would also add a paragraph explaining in more details the construction of the crop rotations and their interactions.

Details about the management are missing. How are the dates for the different agricultural management (sowing, tillage, fertilizer application, etc.) options are derived? Are they constant or dynamic? It is referred to farmers knowledge, but it would be good to get an idea about temporal and spatial variability.

The crop rotation was static for all polygons. The system has been tested before for single polygons. Timings and management details have been derived from local farmer advisers.

Why is the yield not evaluated annually?

We aimed with the study to address the presentation of the full N budget of a cropping system. The reporting and evaluation of the yield is only to conclude on the models performance against the only available validation data. Again, the uncertainty in the reported data is high as e.g. data for food corn and silage corn is aggregated into one category…