Interactive comment on “ Uncertainties in model predictions of nitrogen fluxes from agro-ecosystems in Europe ” by J . Kros

We thank the reviewer for the appraisal and suggestions, which have helped us to improve the manuscript. We have carefully tried to address the issues raised and to revise the paper accordingly. In the revision we have been able to incorporate nearly all the suggestions of the referees as explained in our responses to each reviewer. Below the essence of the questions and suggestions of the reviewer (RE) are given along with our author (AU) replies


Introduction
The nitrogen cycle is of fundamental importance in ecosystem functioning, global change and human health issues. Nitrogen provides a key control of the global carbon cycle through effects on primary production and decomposition.
It is a major determinant of terrestrial and aquatic biodiversity (Dise et al., 2011). It influences particle formation and other chemical production processes in the atmosphere and has major impacts on greenhouse gas (GHG) fluxes via nitrous oxide (N 2 O) and indirectly carbon dioxide (CO 2 ) and methane (CH 4 ) (Butterbach-Bahl et al., 2011;De Vries et al., 2011a).
Several estimates have been made of European-wide land nitrogen (N) budgets, including emissions of ammonia (NH 3 ), nitrogen oxides (NO x ) nitrous oxide (N 2 O) and di-nitrogen (N 2 ), and the leaching and runoff of N (mainly nitrate, NO 3 − ) to ground-and surface waters (Stohl et al., 1996;Bouwman et al., 1997;Mosier et al., 1998;Freibauer, 2003;Van Drecht et al., 2003). These estimates were based on coarse or generic data, because much of the information required for a detailed assessment of N emissions from agriculture is not included in current European databases (Leip, 2004). De  disclosed these detailed data, which are used in the model INTEGRATOR for Europeanwide N and GHG emissions estimates. The model was developed within the framework of the European Integrated Project NitroEurope (Sutton et al., 2007) to model the response of reactive nitrogen (Nr) and GHG emissions due to land-use changes  and the effect of mitigation measures .
Until now, a systematic uncertainty assessment of all N fluxes from terrestrial systems at large regional scales is lacking almost completely and uncertainties reported in the literature vary widely. One of the few examples is from De , who performed an uncertainty analysis of all major N flows in the Netherlands. This research, however, neglected the uncertainty due to spatial model inputs such as soil type and land use, which also influence the uncertainty in GHG predictions (Mosier et al., 1998;Pihlatie et al., 2004) and spatial correlations between these uncertain model inputs. Most research has been focussed on uncertainties in N 2 O emissions, which have been quantified at field scale (e.g. Lehuger et al., 2009), at landscape scale (e.g. Nol et al., 2010) and at national scale (e.g. Del Grosso et al., 2010). Nol et al. (2010) performed an uncertainty propagation analysis for a Dutch fen meadow landscape using the INITIATOR model, while including uncertain spatial model inputs. They used a Monte Carlo (MC) analysis combined with a novel method for estimating and simulating continuous-numerical and categorical input variables, handling spatial and crosscorrelations. Del Grosso et al. (2010) performed a nationwide uncertainty analysis, using the DAYCENT model combined with an empirical approach, to quantify uncertainties in soil N 2 O emissions from croplands in the USA. Uncertainty assessments have also been done for total atmospheric emissions of NH 3 , N 2 O and NO x for various countries, such as the Netherlands (Van Gijlswijk et al., 2004).
Apart from uncertainty assessments that focus on the propagation of model input uncertainty to model outputs, insight in the uncertainty of N fluxes can also be derived from model inter-comparisons. Results of INTEGRATOR, for example, have been thoroughly compared with the results of other European-wide N budget models, including NH 3 , N 2 O and NO x emissions and the sum of N leaching and N runoff to ground-and surface water . Furthermore, Leip et al. (2011b) provided an overview of European-wide N 2 O estimates based on different empirical and process-based model applications.
The aim of this paper is to analyse how uncertainties in model inputs and parameters propagate to model outputs, focussing on uncertainties in continuous model inputs (e.g. fertilizer use) and model parameters (e.g. excretion factors and emission fractions). We present an MC propagation analysis with INTEGRATOR at the European scale, while taking spatial correlation into account. The considered outputs are N 2 O, NO x and NH 3 emission, N leaching to groundwater and N runoff to surface water for the year 2000. Uncertainties in climate, land cover and soil type were not included. The study focussed at (i) uncertainty quantification (UQ) of N fluxes for the entire EU27 and for individual member states, (ii) the identification of the main sources of uncertainty for all N output variables (uncertainty analysis, UA), and (iii) the quantification of the robustness of the uncertainty assessment (robustness analysis, RA).

Modelling N fluxes in agriculture at the European scale
INTEGRATOR uses (i) relatively simple and transparent model calculations based on existing model approaches, and (ii) high-resolution spatially explicit input data. INTE-GRATOR includes various sub-models for the prediction of N (NH 3 , NO x , N 2 O and N 2 ) emissions and N leaching from (i) housing and manure storage systems and agricultural soils, i.e. an adapted version of the MITERRA-Europe model (Velthof et al., 2009;Lesschen et al., 2011), (ii) non-agricultural terrestrial systems, and (iii) an emissiondeposition matrix for NH 3 and NO x , based on the EMEP model (Simpson et al., 2006) to assess interactions between agricultural and non-agricultural land. INTEGRA-TOR ( Fig. 1) calculates the total manure production for each FSSNUTS region, i.e. Farm Structure Survey (FSS) regions which are either at NUTS (Nomenclature of Territorial Units for Statistics) 2 or 3 level (Leip et al., 2008). The manure production is calculated from the animal numbers, available at FSSNUTS level, and the N excretion per animal category, available at country level. A division is made between excretion of animals in housing systems and grazing animals in pastures, based on data at country level. All manure produced in housing and manure storage systems within an FSSNUTS region, corrected for nutrient losses (gaseous and leaching) in these storage systems, is applied to agricultural land in the same region (no manure transport between regions). Weighing factors for grassland, fodder crops and other crops are used to distribute animal manure, which are country dependent (Velthof et al., 2009). Fertilizer N application is based on national fertilizer consumption rates for the year 2000 (FAO, 2010). For each country, the mineral fertilizer is distributed over crops using weighing factors that are based on the N uptake of the crop (sum of N in harvested products and N in crop residues). N uptake is calculated as the product of yield and N content, where the N content depends on the current N input. The yields of arable crops for each country were derived from FAO statistics on a country basis (FAO, 2010). The N contents of harvested crop products and the amount of crop residues were derived from the literature (see Velthof et al., 2009).
The emission of gaseous N compounds (NH 3 , N 2 O, NO and N 2 ) accounted for in the model includes emission (i) from faeces and urine during storage in housing and manure storage systems, (ii) by grazing animals, (iii) after application of manure and fertilizers to agricultural land and (iv) due to atmospheric deposition, nitrogen fixation and crop residue input (not included for NH 3 ). NH 3 , N 2 O and NO x emissions from housing/manure storage and grazing are calculated by multiplying the number of animals and excretion rates per animal category by country-specific emission factors based on the GAINS model (Klimont and Brink, 2004). NH 3 emissions due to fertilizer and animal manure application are calculated by multiplying the application rates of different animal manures and fertilizers by animal and country-specific (related to application techniques) emission factors based on the GAINS model (Klimont and Brink, 2004). Contrary to MITERRA-Europe, which uses generic N 2 O emission factors for fertilizer and animal manure application, INTEGRATOR uses soil emission factors that vary with the N source (three fertilizer types, seven manure types, three crop residue types, mineralized soil organic N, biological N fixation and atmospheric deposition), manure application technique, soil type, land use and precipitation (Lesschen et al., 2011). NO emissions are linearly related to N 2 O emissions.
Losses of N from agricultural systems to ground-and surface waters accounted for in the model include (i) leaching from stored manure to groundwater, (ii) surface runoff to surface waters, (iii) subsurface runoff to surface waters, and (iv) leaching to groundwater (Fig. 1). Leaching from stored manure to groundwater is calculated by multiplying the amount of N stored by leaching fractions that depend on the storage system (liquid or solid manure) and the presence of concrete floors. Data on the distribution of the storage systems and the percentage of covered manure storage in countries were derived from GAINS (Klimont and Brink, 2004). Surface runoff is calculated as a fraction of the various N inputs, using runoff fractions that depend on slope, precipitation, land cover, soil type and soil depth (Velthof et al., 2009). The sum of N leaching and subsurface runoff from soils is derived by multiplying the N surplus by leaching fractions, which are determined based on soil texture, land use, precipitation surplus, soil organic carbon content, temperature and rooting depth (Velthof et al., 2009). The remaining fraction is assumed to be denitrification to N 2 . The N surplus in soils is calculated from the total N input, corrected for gaseous N losses and surface runoff losses directly following application, and the N removal via harvested crops. The division of total N losses to water over leaching and subsurface runoff is based on the water fluxes in both directions, which is based on Keuskamp et al. (2012). Indirect N 2 O emissions from N leached to ground-and surface waters are calculated using IPCC (2006) emission factors.
Three types of grasslands (intensively and extensively managed grasslands, and rough grazing), 30 crop types and eight animal categories (dairy cows, other cows, pigs, laying hens, other poultry, horses, sheep and goats and fur animals) are distinguished. Animal numbers for the year 2000 per NUTS2 were taken from GAINS (Klaassen et al., 2004), and (area weighted) downscaled to the FSSNUTS regions. Data on housing systems and grazing periods were also derived from GAINS, based on national data provided by national experts. Emissions are calculated for spatial units that consist of clusters of 1 km 2 grid cells (NCU, NitroEurope Computational Units), characterized by similar environmental and/or agronomic conditions. For the EU27 (excluding Malta and Cyprus), 35 101 NCUs were distinguished and 744 FSSNUTS regions. To avoid confusion with the EU25 established in 2004, i.e. EU27 without Bulgaria and Romania, we used "EU27" to depict the analysed member states in this study.

Overall approach
Model output uncertainty is determined by three categories of uncertainty sources: (1) model input uncertainty, (2) model structure uncertainty, and (3) model solution uncertainty. In this paper we solely focussed on model input uncertainty. We assumed that model solution uncertainty, which refers to errors caused by rounding, numerical evaluation of integrals, suboptimal optimization solutions, etc., has a marginal contribution to the output uncertainty and can therefore be ignored. Uncertainties due to model structure are not easy to quantify. A possible approach is to compare results of INTEGRATOR with results from other models, which has been done by De .
In this study we performed both an uncertainty quantification (UQ) and uncertainty analysis (UA) using a Monte Carlo (MC) simulation approach. Attractive properties of the MC method are easy implementation, general applicability and resulting in an entire probability distribution of the model output (Heuvelink, 1998). The UQ analyses how uncertainties in model inputs and model parameters propagate to the model output, while UA quantifies the contribution of individual sources of uncertainty to the output uncertainty. The application of UQ and UA involved the following steps: 1. Select model input parameters to be included in the uncertainty propagation analysis.
2. Parameterize the probability distribution of the uncertain inputs.
3. Generate realizations of the uncertain inputs by random drawing from their probability distributions.
4. Submit realizations to INTEGRATOR and perform MC runs.
5. Calculate summary measures of the outputs of the MC runs (e.g. mean, standard deviation, confidence limits).
For UQ we considered all model inputs parameters as uncertain, ranging from activity data to model constants. For UA we selected groups of model parameters for which we estimated the uncertainty contribution. INTEGRATOR was applied while using the NCU as the spatial support, implying that spatial variability within an NCU was not taken into account. The size of the NCUs was highly variable with a mean area of 163 km 2 and standard deviation of 557 km 2 . Larger NCUs occur mainly in flat areas with uniform soil types, whereas small NCUs occur in mountainous areas with variable soils.

Assessment of the uncertainty of the INTEGRATOR model inputs
The UQ focussed on analysis of propagation of uncertainty in (1) initial values (i.e. values of state variables at the start of the simulation), (2) model parameters and (3) environmental constants and variables. In the sequel these three groups are referred to as "model inputs". Uncertainties in categorical data such as land use and soil maps were not considered. The uncertainty in model input data was limited to model inputs affecting (i) N inputs to the system, i.e. N fixation, N deposition, N manure input and N fertilizer, and (ii) N fluxes in and from the ecosystems. For each model input included in the UQ we defined the following statistical properties: 1. Uncertainty in terms of coefficient of variation (CV) or standard deviation (SD), distribution type (normal or lognormal), minimum (min) and maximum (max) at NCU level.
2. Cross correlation for certain pairs of model inputs i and j (ρ cc (i,j )).
3. Spatial correlation coefficients between locations within the different spatial scale levels.
In assessing the uncertainty of model inputs and their spatial correlation, we distinguished four spatial scale levels in order of increasing size: (i) 35 101 NCUs characterized by similar environmental and agronomic conditions, (ii) 744 FSSNUTS regions, (iii) 25 countries, and (iv) EU27 as a whole. An example of these spatial levels is presented for the Netherlands and surroundings in Fig. 2. For model inputs with assumed lognormal distributions; the spatial and cross-correlations of these variables were defined at the log-transformed scale. We used a lognormal distribution when data showed that the variable has a skew distribution. The tails of the distribution can have a marked influence on model outcomes and it is important to take this into account. Details of each aspect are further described below.

Values used for coefficient of variation and standard deviation
To characterize the uncertainty in model inputs, we used the CV for normally distributed model inputs, assuming that the SD is proportional to the average value of a model input. For lognormal model inputs we used the conventional approach and defined the uncertainty by the SD at the log-transformed scale, being approximately equal to the CV on the original scale (Limpert et al., 2001).
Because we have little quantitative information on the uncertainty of the model inputs, we decided to use only three levels of CVs, i.e. low (CV = 0.10), moderate (CV = 0.25) and high (CV = 0.50), with the following assignments: -Low -used for model inputs based on accurate quality statistics for agronomic data (i.e. animal numbers, national N fertilizer inputs and N input data) and also for the crop uptake efficiency fractions.
-High -used for guestimated model inputs, such as N 2 O emission fractions from solid manure systems and N fixation for arable land and grassland.
-Moderate -used for all other model inputs.
The normal distribution was used as a standard, except for (i) all N 2 O and NO emission fractions, and (ii) the ratio between NO x and N 2 O emission fractions. These model inputs, with a strongly positive skew distribution, were considered lognormally distributed.
Since the information on the assigned CVs was largely based on expert judgement, we applied a robustness analysis by using three uncertainty scenarios (Optimistic (O), Reference (R) and Pessimistic (P)). For the CVs, explicit values were assigned for each scenario (see Table 1).
In order to prevent physically unrealistic values, we set (i) a minimum of 0 and a maximum of 1 in case of fractions, and (ii) in other cases a physical minimum (generally 0, sometimes minus infinity) and a physical maximum (generally infinity).

Cross-correlations
Some model inputs are correlated with other model inputs, which affects the uncertainty estimation; therefore crosscorrelations need to be considered. The considered crosscorrelations are given in Table 2. The selected pairs of correlated inputs and assigned correlations were determined by expert judgement. It is obvious that the N excretion rate of cattle is positively correlated with the N content of grass intensive. The negative correlation between Crop yield and the Maximum N content in the harvested crops prevents the occurrence of unrealistic values of crop removal (in terms of N). Although processes related to N 2 O emission and NH 3 emission are different, we also included a correlation between the NH 3 emission fraction from housing systems and those for N 2 O emission fraction from housing systems. This was done since "open" stables or non-covered manure storages leading to a high NH 3 emission also result in a higher N 2 O emission. Since NO and N 2 O emission from housing systems are highly interconnected processes, we assigned a high correlation between the emission fractions for these N forms.

Spatial correlations
In addition to cross-correlations, spatial correlation of uncertain inputs was considered as well. Uncertainty about spatially distributed inputs tends to be positively spatially correlated, and this influences the degree to which uncertainties cancel out by spatial aggregation (see e.g. Heuvelink and Pebesma, 1999). The common geostatistical procedure to include spatial correlations in uncertain inputs is to define variograms and determine for each model input the sill, nugget, range and shape of the variogram. Cross-variograms are also needed for inputs that are cross-correlated. Since hardly any data are available to derive these variograms and cross-variograms, we included spatial correlation using a pragmatic approach as suggested by Lesschen et al. (2007). They incorporated the degree of spatial correlation as an effect on the variance of the aggregated results. The methodology by Lesschen et al. (2007) is also attractive because (i) it can easily be applied in a situation with variable support, whereas the geostatistical procedure assumes a constant support, and (ii) it can take explicit differences between countries into account, e.g. due to legislation and/or country-specific management, whereas the geostatistical procedure cannot cope easily with administrative boundaries. The degree of spatial correlation was defined as the correlation between inputs in (i) different NCUs within the same FSSNUTS region (ρ NCU ), (ii) different FSSNUTS regions within the same country (ρ NUTS ), and (iii) different countries within EU27 (ρ COUNTRY ). The uncertainty of inputs at NCU level was assumed to be constant and independent of the size of the NCU. The uncertainty of most inputs, however, was only available at either FSSNUTS level, e.g. animal numbers, or at national level, e.g. fertilizer amounts and model inputs on N excretion, distribution and emission. For model inputs whose uncertainty was defined at FSSNUTS or country level, the value of ρ NCU was set to 1 (perfect correlation). We also assumed that there was no spatial variability within NCUs, thus correlations within an NCU were always 1.
To complete the definition of uncertainties in spatial inputs with probability distributions, the next step was to specify the spatial cross-correlation coefficients between a model input at some location with another model input at another location. These were calculated as where i and j refer to model inputs i and j and where ρ cc (i,j ) is the correlation coefficient between model inputs i and j for the same plot. We limited the number of values for the spatial correlation coefficients ρ to five classes: -Perfect -for model inputs that are not linked to the NCU but to a higher aggregation level.
-High -in case of a strong spatial correlation.
-Moderate -in case of a moderate spatial correlation.
-Low -in case of a weak spatial correlation.
-None -in case of no spatial correlation.
Furthermore, as with the CV, we used three scenarios (Optimistic (O), Reference (R) and Pessimistic (P)) to investigate the robustness of the assigned correlations. The values of the  spatial correlation coefficients in all these cases are given in Table 3.

Generation of multiple model inputs
We generated sets of multiple model inputs randomly drawn from the input distribution functions, using the rmultnorm function from the MSBVAR library in R (http://www. r-project.org/) and linking each model input to the NCU, NUTS, Country or EU27 level. Spatial correlation coefficient matrices were built by linking all model inputs to each of the 35 101 NCUs, while incorporating the ρ NCU , ρ NUTS and ρ COUNTRY coefficients. At the same time the crosscorrelations were taken into account. Based on preliminary testing, we decided to perform 1000 MC runs. With this relatively high number of MC runs, the predefined distributions and their correlation structure were adequately represented. The analysis was done for the three robustness scenarios, so in total three sets of 1000 MC runs were generated.

Quantification of the uncertainty contribution
The INTEGRATOR model includes a large number of model inputs, which hampers the quantification of the uncertainty contribution of individual parameters. Therefore, we decided to quantify the contribution of model inputs to the uncertainty in model outputs for groups of model inputs affecting certain model outputs, such as N excretion, N emission, N uptake, etc. We grouped the INTEGRATOR inputs into six groups for which we analysed the uncertainty separately (see column "Group" in Table 4). For each individual model input the uncertainty contribution to the total output uncertainty was estimated. This was accomplished by a MC simulation considering only the model input group uncertain for which the contribution was to be estimated. The other six groups were considered certain by using their default values as stored in the INTEGRATOR database. Hence, the model was run again in six separate rounds, and each time only one of the six groups was made uncertain. The uncertainty contribution was based on a comparison of the resulting output variances of each model input group. The analysed model inputs (51 in total) and their statistical properties and spatial levels are given in Table 4.

Uncertainty at EU27 and country level
Results at EU27 level show larger uncertainties for N leaching to groundwater and N runoff to surface water (relative errors of ∼ 19 %), and smaller uncertainties for the emission of N 2 O, NO x and NH 3 (relative errors of ∼ 12 %) ( Fig. 3 and Table 5).
The larger uncertainties in N runoff to surface water and N leaching to groundwater are caused by the larger number of uncertain model inputs affecting excretion, emission, nitrogen input, nitrogen uptake and leaching and runoff, whereas the uncertainties in N 2 O, NO x and NH 3 are affected by fewer uncertain model inputs, i.e. excretion and emission (see Table 4).
The uncertainties can vary considerably among countries (Table 6). Some countries have large uncertainties for all outputs, e.g. Austria, whereas other countries have relatively high uncertainties for only one output, e.g. Denmark for N leaching to groundwater. As shown in Table 5, the uncertainties in the average emissions and leaching fluxes at European level are noticeably smaller than those at country level. The average relative uncertainty for N 2 O at country level is 16 % (Table 6), whereas at EU27 level it decreases to 12 % (Table 5).

Spatial distribution of uncertainty at NCU level
Figures 4 and 5 depict the spatial variability of the areaweighted average N (N 2 O, NH 3 and NO x ) emissions (Fig. 4) and N losses to groundwater and to surface water (Fig. 5) within the EU27 at NCU level for the year 2000. The figures show the mean values and the relative uncertainty expressed as a CV. The results show that N fluxes are highly correlated with agricultural intensity. High fluxes occur in areas with a high intensity of animal husbandry, such as the Po Valley, northern Germany, the Netherlands and Bretagne, and hilly regions such as the Alps in southern Germany and northern Spain (see e.g. Leip et al., 2011a). Compared to other outputs, NH 3 emissions show a smaller spatial variation. NH 3 emissions are related to livestock manure production and manure application, which are provided at FSS-NUTS level, whereas N 2 O and NO x fluxes and N leaching and runoff fluxes also are largely influenced by the underlying environmental conditions such as soil type and hydrological conditions, which are linked to the NCU level.
The variation in the uncertainty of the N 2 O and NO x emission among the NCUs is relatively small, and varies mostly between 20 and 40 %. The uncertainty at country level is lower, generally below 20 % for both N 2 O and NO x (Table 6).
The spatial variation of the uncertainty of NH 3 emission at NCU level is somewhat larger and apparently distributed into two classes, either less than 20 % or between 20 and 40 %. It seems that areas with a higher uncertainty are areas where most of the NH 3 emission is caused by grazing and manure application, such as in Spain and Ireland.
Contrary to the gaseous emissions (Fig. 4), N runoff to surface water and especially N leaching to groundwater show a relatively large spatial variation (Fig. 5). In large parts of Central and Northern Europe, but also in parts of the UK, Ireland, France, Spain and Italy the uncertainty at NCU level of N leaching to groundwater is even larger than 100 %. For N runoff to surface water there are fewer areas with an uncertainty of more than 100 %. Large uncertainties in N leaching to groundwater generally occur in countries with a relatively large area of sandy soils, for which the uncertainty is large compared to clay and peat soils, as sandy soils generally receive higher N inputs and are more susceptible to leaching (not shown). Moreover, the large variation in the uncertainty for leaching and runoff is mostly due to the leaching fraction and runoff fractions, which are related to soil type, texture and organic matter.
Uncertainties increase when going from EU27 level to country level to NCU level, as shown in Table 7. This table shows that the uncertainties (expressed as CV) for the outputs at NCU level are ∼ 3 times greater, and at country level    Table 1. 5 See Table 3. 6 We assigned the uncertainty to each applicable emission factor, related to the different N sources, i.e. deposition, fixation, animal manure, fertilizer, mineralisation and crop residues. 7 In INTEGRATOR, NO x soil emissions were derived as a fraction of the N 2 O soil emissions. The CV of this fraction was based on the dataset from Stehfest and Bouwman (2006). For this parameter we used 0.5 × CV for the O scenario, CV for the R scenario and 1.5 for the P scenario.

Uncertainty contribution of various inputs
Uncertainty in the N 2 O emission is mainly caused by uncertainty in N inputs and crop uptake parameters (each group contribute > 20 %) and to a lesser extend to the uncertainty in housing emission (19 %), excretion (10 %) and soil emission parameters (2 %) (Fig. 6). The same is true for the NO emission, but the uncertainty contribution of the soil emission parameters (5 %) is somewhat larger and the contribution of housing emission parameters is smaller (4 %). For NH 3 emissions, the largest uncertainties originated from the excretion parameters (51 %), followed by the uncertainty in N inputs (20 %). As with the N 2 O and NO emissions, the uncertainty contributions of the housing emission (14 %) and soil emission parameters (2 %) are rather small. The uncertainty in N leaching and N runoff is mainly caused by uncertainty in N inputs (48 % and 57 % resp.), N leaching parameters (26 % and 17 % resp.) and crop uptake (13 % and 16 % resp.).

Robustness analysis
The results of the robustness analysis are presented in Table 8. Results show that the uncertainty in N losses to air (N 2 O emission, NO x emission, NH 3 emission) or water (N leaching to groundwater and N runoff to surface water) for the EU27, expressed as CV, range from 4 to 6 % for the Optimistic scenario, from 12 to 19 % for the Reference scenario and from 19 to 30 % for the Pessimistic scenario. The maximum values of the Optimistic and Pessimistic scenarios range from 12 to 55 %. The uncertainty due to the three robustness scenarios is about 50 % (CV rob , Table 8).  The robustness analysis shows that the mean N losses also depend on the scenario. The mean values increase when shifting from the Optimistic to the Pessimistic scenario. This is due to non-linearities and the lognormal distribution of some of the uncertain model inputs, which causes a shift towards higher values.

Limitations of the study
We have used a methodology allowing inclusion of different degrees of spatial correlation in uncertain model inputs. The parameterization, i.e. both the quantification of the statistical properties and the linkage of model inputs to the spatial levels and groups, is difficult because much of the required information is lacking. Nevertheless, we tried to incorporate most of the uncertainties and their correlations. It is evident that all assumptions made influence the results. In order to compensate for these shortcomings, we included a robustness analysis with respect to the assigned uncertainties. Yet, future research should address the improvement of the parameterization of uncertainty distributions and their (spatial) correlations, but this requires that appropriate data become available. Other aspects that were not included in the uncertainty quantification deserve attention too, such as uncertainty due to model structure, and uncertainty in categorical input data, such as land use maps and soil maps.

Plausibility of the uncertainty quantification
Leaching fluxes to groundwater and runoff fluxes to surface water have the largest relative (CV) and absolute (SD) uncertainties. NH 3 emission has the smallest relative uncertainty (CV), but the largest absolute uncertainty (SD). In general the results show that the further in the calculation chain, and the more uncertain processes are involved, the larger the uncertainty in the outputs. Furthermore, results confirm that uncertainties and spatial variation in model outputs are partly cancelled out due to spatial aggregation.
Regarding N 2 O emissions, earlier studies show comparable uncertainties. Nol et al. (2010) estimated an average N 2 O emission at landscape scale of 20.5 kg N 2 O-N ha −1 yr −1 and a standard deviation of 10.7 kg N 2 O-N ha −1 yr −1 , implying a relative uncertainty of 52 %. This is clearly larger than the uncertainty at the EU27 level from this study (12 %), but also larger than uncertainties from this study at NCU level (26-44 %). The reason for the larger uncertainty might be that Nol et al. (2010) also included uncertainty on detailed management information, which is lacking in our study. Due to the different nature and character of the agricultural sources, among regions and countries, as well as by the limited number and uneven spread of the measurements concerning related emission parameters, the uncertainty in local scale estimates of N 2 O emissions can be rather large .
At the national scale De  estimated the uncertainty of modelled N fluxes from agriculture for the Netherlands. At national level they quantified a relative uncertainty of 24 % for N 2 O, 16 % for NO and 16 % for the NH 3 emissions, and 47 % for the leaching and runoff to ground-and surface water. Except for NH 3 , these uncertainties are 20 to 50 % larger than found in our study for the Netherlands (see Table 6). Since De  ignored spatial correlation by assuming a perfect correlation between the uncertainties of their spatial units, it is likely that De  overestimated the uncertainty in the national N fluxes.
From a global scale uncertainty analysis, Beusen et al. (2008) calculated a range of 27-38 (with a mean of 32) Tg yr −1 for the global NH 3 -N emission from agricultural systems, i.e. a relative uncertainty of 11 %, which is slightly lower than the 16 % from this study. In addition, Van Gijlswijk et al. (2004) performed an uncertainty assessment for the total atmospheric emissions of NH 3 using a Tier 1  an uncertainty range of about −35 % to +50 % or a relative uncertainty of 42 % at the national level. The relative uncertainty tended to be larger at the regional level, particularly in regions with low emissions. The national level uncertainty for the USA is clearly larger than the 17 % from this study for the EU27 level. However, the study from Del Grosso et al. (2010) also included model structure uncertainty, which contributed 83 % to the total uncertainty. The remaining 17 % was due to uncertainty in model input data, implying a relative uncertainty of 7 % (0.17 × 42 %). This is smaller than the 12 % from our study. The difference might be caused by the fact that the Del Grosso study is confined to soil emission, whereas this study also included housing emissions. Moreover, the model structure uncertainty of Del Grosso et al. (2010) also included parameter uncertainty, which was considered part of the model input uncertainty in this study. For the EU25, Schulze et al. (2009) estimated an uncertainty in total N 2 O emission from agriculture of 50 %, which was based on the IPCC guidelines 2006 (IPCC, 2006) while assuming a perfect spatial correlation. Based on the results of this study and other studies discussed here, it is likely that Schulze et al. (2009) overestimated the uncertainty in the N 2 O emission. Moreover, Leip (2010) thoroughly discussed the uncertainty of N 2 O emission estimates in GHG inventories and demonstrated the importance of correlation, if uncertainties are combined for the whole of Europe. He concluded that the uncertainty estimates tend to be underestimated for the activity data and overestimated for the emission factors. Berdanier and Conant (2012) showed that improved regional emission factors will reduce uncertainty by up to 65 % for estimates of regional and global N 2 O emissions from crop land.
The plausibility of the INTEGRATOR results was examined by De Vries et al. (2011b) by comparing the computed mean values against estimates by three other models with different complexity and data requirements. They found that the estimated overall variation at EU27 level is small for the emissions of ammonia and nitrous oxide, but large for N leaching and runoff. At smaller spatial level, however, large differences in N output fluxes were found.
The analysis presented in our paper was limited to input uncertainty and ignored uncertainty due to model structure. However, model structural uncertainties can be large because these include the many assumptions and simplifications that any modelling exercise makes (see e.g. Del Grosso et al., 2010;Refsgaard et al., 2006;van der Sluijs, 2007). Possibilities for better insight in model structure uncertainty are to put more emphasis on the inter-comparison of results from independently developed models  and by using a validation with independent observations.

Contribution of different sources to the uncertainty
The uncertainty of most of the considered outputs is mainly influenced by the uncertainty in model inputs that influence the N soil input such as the allocation of animal manure and fertilizer use. Notable is the relatively large contribution due to N uptake (10-30 %), whereas the contribution of emission factors (for housing and soil) both for N 2 O and NH 3 is rather small (0-20 %). Apparently, their contribution is overruled by other uncertainties in the N cycle. This implies that more emphasis on collecting reliable N uptake data can be more helpful to reduce the uncertainty in the considered N fluxes for EU27, rather than putting more effort into improving emission factors. However, this low contribution of the emission factors is partly due to the choices made on the spatial correlation (cf. Table 4), which might be disputed. A robustness analysis of the uncertainty analysis could be recommended as well. For N leaching, the uncertainty contribution of the leaching parameters is larger compared to the gaseous N emissions parameters. Improving N leaching and runoff fractions should therefore receive more attention.
The uncertainty contribution was based on a comparison of the resulting output variances, which has some disadvantages, as the addressed uncertainty is only based on the inputs that are considered uncertain. The other inputs have fixed reference values that are not necessarily true values. In addition, the correlations between the uncertain and the fixed inputs are ignored. Although these disadvantages cannot be completely avoided, there are techniques developed to tackle the problem more thoroughly (see e.g. Jansen et al., 2000;Saltelli et al., 2004), but these techniques are computationally more intensive.

Conclusions
This study shows large variation in uncertainties and uncertainty contributions of the different model input sources across the various model outputs. Furthermore, the uncertainty in nitrogen fluxes differs per country.
The most prominent result of this study is that the uncertainty at the EU27 level is smaller than 30 % for all considered model outputs, which is rather low. The results indicate that the uncertainty in all N fluxes at EU27 level is most likely in the range of 5-30 %, including all three scenarios of the robustness analysis,and increases in the direction NH 3 < NO x < N 2 O< N leaching and runoff. Using the reference scenario, the uncertainty in N fluxes varies from 11 to 19 %. The relative uncertainty in model outputs is 1.3-2 times larger for the pessimistic scenario and ranges from 16 to 30 %, while the uncertainty is 2-3.3 times smaller for the optimistic scenario and ranges from 4 to 6 % only. The spatial variation in the uncertainty in the N fluxes can be large, in particular for the leaching fluxes to groundwater (N le gw ).
Uncertainties in N fluxes at EU27 level were smaller at larger spatial scales than at smaller spatial scales. This is in line with expectations because uncertainties in model outputs are partly cancelled out at higher scale levels due to spatial aggregation. Although an objective comparison is not straightforward, the uncertainties in N fluxes at EU27 level are generally smaller compared to results from other studies. Presumably most regional uncertainty assessment studies ignore spatial correlation and assume that inputs are constant over space, which leads to an overestimation of uncertainty at coarse spatial scales.
Activity data, such as N fertilizer, N manure and N uptake, are the main sources of uncertainty. It is advised to put more effort into obtaining reliable information about those model inputs, rather than putting more effort into improving the accuracy of emission factors.
As far as we know, this research is the first attempt to quantify the uncertainty in N fluxes at the EU27 level, while including a correlation between spatially variable model inputs, which is crucial in regional scale assessment. Unfortunately, for many of the spatially variable inputs very little information on their uncertainties was available. Therefore, the uncertainty of most of the model inputs was based on guestimates rather than real data. Further effort should be put into the gathering of spatial variable input data and quantification of associated uncertainties.