# Long-term annual soil nitrogen surplus across Europe (1850–2019)

Masooma Batool, Fanny J. Sarrazin, Sabine Attinger, Nandita B. Basu, Kimberly J. Van Meter, Rohini Kumar

► To cite this version:

Masooma Batool, Fanny J. Sarrazin, Sabine Attinger, Nandita B. Basu, Kimberly J. Van Meter, et al.. Long-term annual soil nitrogen surplus across Europe (1850–2019). Scientific Data, 2022, 9 (1), pp.612. 10.1038/s41597-022-01693-9. hal-04760886

HAL Id: hal-04760886

https://hal.science/hal-04760886v1

Submitted on 31 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

DATA DESCRIPTOR

<figure><img src="image_2.png" /></figure>

Check for updates

# Long-term annual soil nitrogen surplus across Europe (1850–2019)

Masooma Batool<fnref n="1" />✉, Fanny J. Sarrazin<fnref n="1" />, Sabine Attinger<fnref n="1" />,², Nandita B. Basu³,⁴,⁵, Kimberly Van Meter⁶,⁷ & Rohini Kumar<fnref n="1" />✉

Worldwide surface waters suffer from the presence of nitrogen (N) compounds causing eutrophication and deterioration of the water quality. Despite many Europe-wide legislation's, we still observe high N levels across many water bodies in Europe. Information on long-term annual soil N surplus is needed to better understand these N levels and inform future management strategies. Here, we reconstructed and analysed the annual long-term N surplus for both agricultural and non-agricultural soils across Europe at a 5 arcmin (≈10 km at the equator) spatial resolution for more than a century (1850–2019). The dataset consists of 16 N surplus estimates that account for the uncertainties resulting from input data sources and methodological choices in major components of the N surplus. We documented the consistency and plausibility of our estimates by comparing them with previous studies and discussed about possible avenues for further improvements. Importantly, our dataset offers the flexibility of aggregating the N surplus at any spatial scale of relevance to support water and land management strategies.

## Background & Summary

Worldwide, surface waters suffer from high concentration of nitrogen (N) compounds, mainly due to large usages of agrochemicals (inorganic fertilizers)<fnref n="1" />,² and deposition of atmospheric N originating from fossil fuel combustion³⁻⁵. This has led to the deterioration of the surface and subsurface water quality, causing blooms of harmful algae and threatening aquatic ecosystems and human health when N compounds are present in drinking water⁶⁻⁹. To address these N related problems, water quality mitigation measures were implemented, in particular in Europe following the EU Water Framework Directive (Directive2000/60/EC). Despite such measures, we still observe nitrate levels above the regulatory threshold in several water-bodies across Europe¹⁰,¹¹. One of the reasons behind such failure lies in the presence of N legacies, as N levels can depend not only on the current N inputs to the landscape, but also on the past N inputs that have accumulated through time in the soil zone and the groundwater in so-called "legacy stores"¹²,¹³. These legacies are largely responsible for time lags between the implementation of water quality measures and the water quality response¹⁴,¹⁵.

Therefore, information on long-term annual N budgets at relevant scales for water and land management is crucial to better understand N legacies and inform future (land and water) management strategies¹⁶. Existing datasets covering the entire European domain are either based on the N land budget, which is defined as the difference between N inputs (mineral fertilizer, manure, atmospheric deposition, biological N fixation) and N outputs (N from harvested crops and grass removal)¹¹, or the soil surface budget which further excludes gaseous losses during manure management and storage¹⁷. These datasets provide annual N budgets either for limited time periods or at a coarse spatial resolution and they focus on agricultural soils (cropland and/or pasture) only. In particular, the FAO provides global datasets for soil nutrient budgets at country-level for cropland¹⁸ and for different underlying components of the agricultural N budget since 1961. More recently, Zhang et al.⁵ provides estimates for different components of the N budget for cropland, which were based on 13 existing long-term databases covering 115 countries from 1961 to 2015. Further, N budget datasets were developed for European countries specifically, either with respect to total agricultural lands (Eurostat¹⁹; Leip et al.²⁰) or to croplands

<footnote n="1">¹UFZ-Helmholtz Centre for Environmental Research, Department of Computational Hydrosystems, Leipzig, Germany. ²Institute of Environmental Science and Geography, University of Potsdam, Potsdam, Germany. ³Department of Civil and Environmental Engineering, University of Waterloo, Waterloo, ON, Canada. ⁴Department of Earth and Environmental Sciences, University of Waterloo, Waterloo, ON, Canada. ⁵Water Institute, University of Waterloo, Waterloo, ON, Canada. ⁶Department of Geography, The Pennsylvania State University, University Park, PA, USA. ⁷Earth and Environmental Systems Institute, The Pennsylvania State University, University Park, PA, USA. ✉e-mail: masooma.batool@ufz.de; rohini.kumar@ufz.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Workflow for constructing the long-term annual dataset of N surplus during the period 1850–2019. The number in brackets with red colors present different combination of datasets that we used to account for the uncertainties, resulting in 16 N surplus estimates.</figcaption></figure>

only (Einarsson *et al.*²¹). Moreover, efforts have been expanded to provide N budgets at sub-national levels for individual countries, such as, Denmark²², United Kingdom²³, France²⁴, Poland²⁵ and Germany²⁶. Notably, uncertainty estimation is also currently not a common practice in N budget construction¹⁷, although uncertainties in existing N datasets are understood to be large⁵. Uncertainties in the N surplus budget, for example, can stem from a limited knowledge on fertilizer and manure application rates, distinction of their applications among croplands and pasture areas, specifications of livestock and grassland (N) productions, among others⁵,²⁷,²⁸.

To overcome these challenges, we provide a dataset of annual long-term N soil surface budget (termed hereafter as “N surplus”) from both agricultural (cropland and pasture) and non-agricultural soils across Europe at a spatial resolution of 5 arcmin (1/12°; ≈10 km at the equator) over more than a century (1850–2019), while explicitly accounting for uncertainties resulting from input data sources and methodological choices in major components of the N surplus (i.e., fertilizer, manure and N removal rates). This newly developed dataset combines information at different levels (in particular country and grid level) to construct the different N surplus components. We assessed the importance of having a long-term dataset of N surplus, as the magnitude of the N surplus has changed tremendously over the last 100 years. Moreover, we show the relevance of considering N surplus over non-agricultural areas, since these areas account for 30% of the total N surplus on average across the study area of Europe. Our gridded database also offers the flexibility of aggregating N surplus at any spatial scale of relevance for the design of land and management strategies. Importantly, we demonstrate the plausibility and consistency of our N surplus estimates by comparing them with existing N budgets datasets¹¹,²⁰,²²–²⁶,²⁹,³⁰. We then discuss about possible avenues towards a more comprehensive uncertainty characterization in N surplus. In this respect, our reconstruction methodology, described in detail below, paves the way for the exploration of further alternative assumptions in N surplus.

## Methods

This section details our methodology for reconstructing the long-term annual time series of the land use and the individual components of the N surplus at gridded level (5 arcmin) over the time period 1850–2019 (see detailed workflow in Fig. 1). We compiled and harmonized a range of databases available for different time periods (i.e., for the periods 1850–1960 and 1961–2019, and for the year 2000), at different frequencies (i.e., yearly, decadal, snapshots) and at different spatial resolutions (i.e., global trends, country-level values, gridded values) following the harmonization methodologies from previous studies¹,¹⁶,³¹,³². The data sources used in this study are listed in Table 1. Importantly, we accounted for the uncertainties resulting from the differences in underlying data sources and methodological choices in major components of the N surplus. To this end, we constructed 16 gridded time-series of N surplus estimates by combining two estimates for fertilizer, four estimates for animal manure, and two estimates for the N removal from pastures.
<table><thead><tr><th rowspan="2">Datasets</th><th rowspan="2">Variables</th><th colspan="3">Spatial feature</th><th colspan="2">Temporal feature</th><th rowspan="2">Ref</th></tr><tr><th>Extent</th><th>Level</th><th>Resolution</th><th>Extent</th><th>Resolution</th></tr></thead><tbody><tr><td>Ramankutty et al.<sup>34</sup></td><td>Land use (cropland and pasture area)</td><td>Global</td><td>Grid</td><td>5' × 5'</td><td>2000</td><td>Snapshot</td><td>34</td></tr><tr><td>Monfreda et al.<sup>35</sup></td><td>Crop specific harvested area and production</td><td>Global</td><td>Grid</td><td>5' × 5'</td><td>2000</td><td>Snapshot</td><td>35</td></tr><tr><td>HYDE</td><td>Land use (cropland and pasture area)<sup>a</sup></td><td>Global</td><td>Grid</td><td>5' × 5'</td><td>1700–2017</td><td>Decadal or Annual<sup>b</sup></td><td>40</td></tr><tr><td>Global Land Cover (GLC)</td><td>Land use (forest, semi-natural vegetation, urban, non-vegetation area)</td><td>Global</td><td>Grid</td><td>5' × 5'</td><td>2000</td><td>Snapshot</td><td>39</td></tr><tr><td>input4MIPS</td><td>N atmospheric deposition</td><td>Global</td><td>Grid</td><td>1.9° × 2.5°</td><td>1850–2014</td><td>Monthly</td><td>53</td></tr><tr><td>Zhang et al.<sup>29</sup></td><td>Animal manure production</td><td>Global</td><td>Grid</td><td>5' × 5'</td><td>1850–2014</td><td>Annual</td><td>29</td></tr><tr><td>Holland et al.<sup>42</sup></td><td>Fertilizer production</td><td>Global</td><td>Country</td><td>—</td><td>1925–1960</td><td>Annual</td><td>42</td></tr><tr><td>FAOSTAT</td><td>Land use (cropland, pasture), fertilizer application, animal manure, crop harvested area and production</td><td>Global</td><td>Country</td><td>—</td><td>1961–2019</td><td>Annual</td><td>38</td></tr><tr><td>IFA</td><td>Fertilizer application amounts for different crop categories and pasture</td><td>Global</td><td>Country<sup>c</sup></td><td>—</td><td>2014/2015</td><td>Snapshot</td><td>40</td></tr><tr><td>Bayliss-Smith and Wanmali (1984) as cited in Our World in Data (OWD)</td><td>Wheat yield</td><td>Global</td><td>Country</td><td>—</td><td>1850–1960</td><td>Snapshots<sup>d</sup></td><td>57,58</td></tr><tr><td>Einarsson et al.<sup>21</sup></td><td>Estimates for fodder crops (area, production), fertilizer, manure</td><td>Region</td><td>Country</td><td>—</td><td>1961–2019</td><td>Annual</td><td>21</td></tr></tbody></table>

**Table 1.** Datasets used in this study for the reconstruction of N surplus. <sup>a</sup>According to the definition used in our study, pasture corresponds to the land under permanent meadow and pasture. Therefore, we used the category “grazing” of the HYDE dataset to derive the pasture areas in our study. <sup>b</sup>Decadal for the period 1700–2000 and annual for the period 2000–2016. <sup>c</sup>IFA covers estimates for 28 countries (considering European Union (EU-28) as a single country). <sup>d</sup>OWD provide wheat yield for the years 1850, 1909, 1934, 1960.

In particular, we adopted country-level data from FAOSTAT (Food and Agriculture Organization Corporate Statistical Database)<sup>33</sup> that provide the longest data record across a range of variables such as mineral fertilizer, animal manure, crop production and crop harvested area with global coverage during 1961–2019. Furthermore, we utilized the recent country-level dataset of Einarsson et al.<sup>21</sup> available for the period 1961–2019 for animal manure, distribution of fertilizer and manure to cropland and pasture, fodder crop specific harvested areas and fodder crop production. FAOSTAT includes data for total agricultural areas based on national statistics, whereas, Einarsson et al.<sup>21</sup> derives estimates for croplands by combining information from FAOSTAT, Eurostat and several national datasets for European countries. We downscaled the country-level estimates using, for most variables, the spatial variability in land use areas (Ramankutty et al.<sup>34</sup> dataset), and crop specific harvested areas and crop production (Monfreda et al.<sup>35</sup> dataset) for the year 2000. We used consistent grid level datasets<sup>34,35</sup> to avoid a mismatch at grid level due to differences in spatial variation from different data sources. Furthermore we also adapted (annual and decadal, grid-level) information provided by the History Database of the Global Environment (HYDE version 3.2)<sup>36,37</sup> to reconstruct the temporal trajectories of different land cover/uses (e.g., croplands, pastures, forest and other natural vegetated areas).

In the sections below, we first define the N surplus to agricultural and non-agricultural soils. Second, we provide details on the methodology for the reconstruction of the land use areas (cropland, pasture, forest, semi-natural vegetation, urban, non-vegetation) and crop specific harvested area (non-fodder and fodder crops). Third, we further explain the steps adopted to reconstruct the N inputs (through fertilizer, manure, atmospheric deposition and biological N fixation) and N outputs (N removal from cropland, pastures and forest).

**N surplus.** We used the concept of N soil surface balance, which is the difference between N inputs and N output with adjustment of volatilization losses during manure storage system<sup>11,20,26</sup>. Here, total N surplus is the sum of N surplus over agriculture (cropland and pasture) and non-agriculture areas (forest, semi-natural vegetation, non-vegetation and urban) as given in Eq. (1) (all variables are in kg ha<sup>-1</sup> yr<sup>-1</sup>):

$$ \text{Surp}_{\text{soil}}(i, y) = \text{Surp}_{\text{agri}}(i, y) + \text{Surp}_{\text{NonAgri}}(i, y) \qquad (1) $$

where i is grid cell; y is year; Surp<sub>soil</sub> is the total N surplus; Surp<sub>agri</sub> is the N surplus over agriculture area; and Surp<sub>NonAgri</sub> is the N surplus over non-agriculture area. The individual component of N surplus are explained in the following sections.

**N surplus over agriculture soils.** N surplus over agricultural soils accounts for the N surplus over croplands (Surp<sub>cr</sub>) and pastures (Surp<sub>past</sub>). N surplus over croplands and pastures is estimated as the difference between the inputs to cropland and pasture (Inp<sub>cr</sub> and Inp<sub>past</sub>, respectively) through mineral fertilizers (FERT<sub>cr</sub> and FERT<sub>past</sub>, respectively), animal manure (MAN<sub>cr</sub> and MAN<sub>past</sub>, respectively), biological N fixation (BNF<sub>cr</sub> and BNF<sub>past</sub>, respectively) and atmospheric deposition (DEP<sub>cr</sub> and DEP<sub>past</sub>, respectively), and the N outputs (Out<sub>cr</sub> and Out<sub>past</sub>, respectively) through N removal via harvested crops (Rem<sub>cr</sub>) and via animal grazing and cutting of grass (Rem<sub>past</sub>), respectively. These can be summarized in Eqs. (2–8) (all variables are in kg ha<sup>-1</sup> yr<sup>-1</sup>) as:
$$ \mathrm{Surp}_{\mathrm{Agri}}(i, y) = \mathrm{Surp}_{\mathrm{cr}}(i, y) + \mathrm{Surp}_{\mathrm{past}}(i, y) \qquad (2) $$

$$ \mathrm{Surp}_{\mathrm{cr}}(i, y) = \mathrm{Inp}_{\mathrm{cr}}(i, y) - \mathrm{Out}_{\mathrm{cr}}(i, y) \qquad (3) $$

$$ \mathrm{Inp}_{\mathrm{cr}}(i, y) = \mathrm{FERT}_{\mathrm{cr}}(i, y) + \mathrm{MAN}_{\mathrm{cr}}(i, y) + \mathrm{DEP}_{\mathrm{cr}}(i, y) + \mathrm{BNF}_{\mathrm{cr}}(i, y) \qquad (4) $$

$$ \mathrm{Out}_{\mathrm{cr}}(i, y) = \mathrm{Rem}_{\mathrm{cr}}(i, y) \qquad (5) $$

$$ \mathrm{Surp}_{\mathrm{past}}(i, y) = \mathrm{Inp}_{\mathrm{past}}(i, y) - \mathrm{Out}_{\mathrm{past}}(i, y) \qquad (6) $$

$$ \mathrm{Inp}_{\mathrm{past}}(i, y) = \mathrm{FERT}_{\mathrm{past}}(i, y) + \mathrm{MAN}_{\mathrm{past}}(i, y) + \mathrm{DEP}_{\mathrm{past}}(i, y) + \mathrm{BNF}_{\mathrm{past}}(i, y) \qquad (7) $$

$$ \mathrm{Out}_{\mathrm{past}}(i, y) = \mathrm{Rem}_{\mathrm{past}}(i, y) \qquad (8) $$

**N surplus over non-agriculture soils.** N surplus over non-agriculture soils (Surp<sub>NonAgri</sub>) includes N surplus over forest (Surp<sub>For</sub>), semi-natural vegetation (Surp<sub>NatVeg</sub>), urban (Surp<sub>Urban</sub>) and rest non-vegetated areas (Surp<sub>NonVeg</sub>). N surplus over forest is derived as the difference between N inputs to forest (Inp<sub>For</sub>) through biological N fixation (BNF<sub>For</sub>) and atmospheric deposition (DEP<sub>For</sub>) and N outputs (Out<sub>For</sub>) via N removal from forest (Rem<sub>For</sub>). N surplus over semi-natural vegetation (Surp<sub>NatVeg</sub>) covers N inputs via atmospheric deposition (DEP<sub>NatVeg</sub>) and biological N fixation (BNF<sub>NatVeg</sub>). N surplus over the urban and rest non-vegetative areas (Surp<sub>Urban</sub> and Surp<sub>NonVeg</sub>, respectively) is derived as N inputs via atmospheric deposition (DEP<sub>Urban</sub> and DEP<sub>NonVeg</sub>, respectively). These can be summarized in Eqs. (9–15) (all variables are in kg ha<sup>-1</sup> yr<sup>-1</sup>) as:

$$ \mathrm{Surp}_{\mathrm{NonAgri}}(i, y) = \mathrm{Surp}_{\mathrm{For}}(i, y) + \mathrm{Surp}_{\mathrm{NatVeg}}(i, y) + \mathrm{Surp}_{\mathrm{NonVeg}}(i, y) + \mathrm{Surp}_{\mathrm{Urban}}(i, y) \qquad (9) $$

$$ \mathrm{Surp}_{\mathrm{For}}(i, y) = \mathrm{Inp}_{\mathrm{For}}(i, y) - \mathrm{Out}_{\mathrm{For}}(i, y) \qquad (10) $$

$$ \mathrm{Inp}_{\mathrm{For}}(i, y) = \mathrm{DEP}_{\mathrm{For}}(i, y) + \mathrm{BNF}_{\mathrm{For}}(i, y) \qquad (11) $$

$$ \mathrm{Out}_{\mathrm{For}}(i, y) = \mathrm{Rem}_{\mathrm{For}}(i, y) \qquad (12) $$

$$ \mathrm{Surp}_{\mathrm{NatVeg}}(i, y) = \mathrm{DEP}_{\mathrm{NatVeg}}(i, y) + \mathrm{BNF}_{\mathrm{NatVeg}}(i, y) \qquad (13) $$

$$ \mathrm{Surp}_{\mathrm{NonVeg}}(i, y) = \mathrm{DEP}_{\mathrm{NonVeg}}(i, y) \qquad (14) $$

$$ \mathrm{Surp}_{\mathrm{Urban}}(i, y) = \mathrm{DEP}_{\mathrm{Urban}}(i, y) \qquad (15) $$

**Land use categories.** We compiled a range of databases to estimate annual trajectories of coverage areas for different land types including the harvested areas for 17 non-fodder crops and 6 fodder crop categories. These estimates were used in various steps to reconstruct the N surplus, for example to derive crop specific fertilizer application rates, distribution of mineral fertilizer and animal manure to cropland and pastures. In below, we provide details on the derivation of these areas from available databases.

**Cropland, pasture and agriculture area.** Cropland is defined as land used for cultivation of crops including arable crops and land under permanent crops<sup>11,34,38</sup>. Pasture area in our estimates is the land under permanent meadow and pasture<sup>38</sup>. First, to obtain the spatial variability at gridded level, we used the dataset from Ramankutty et al.<sup>34</sup> which provides the global spatial distribution of cropland and pasture area (denoted as A<sub>Ramankutty-crop</sub> (ha) and A<sub>Ramankutty-past</sub> (ha), respectively) at a spatial resolution of 5 arcmin for the year around 2000. The land use data in Ramankutty et al.<sup>34</sup> is derived based on agriculture inventories and land cover data from satellite imageries. To derive the temporal variability in the gridded values of cropland and pasture area, we adopted the land use dataset from the History Database of the Global Environment (HYDE version 3.2)<sup>36</sup>. HYDE provides global cropland and pasture area as decadal values for the period 1700–2000 and annual for the period 2000–2017 at a spatial resolution of 5 arcmin. We generated the annual time series of the cropland and pasture areas from HYDE (denoted as A<sub>HYDE-cr</sub> (ha) and A<sub>HYDE-past</sub> (ha), respectively) for the period 1850–2000 using linear interpolation for every two decadal estimates. For the year 2018 and 2019, we used the same value as of the year 2017. We derived the temporal variability of cropland and pasture area from the HYDE dataset for
each year during 1850–2019, by referencing it to the year 2000 ($R_{\text{HYDE-cr}}$ (-) and $R_{\text{HYDE-past}}$ (-), respectively), as given in Eqs. (16) and (17):

$$R_{\text{HYDE-cr}}(i, y_{1850-2019}) = \frac{A_{\text{HYDE-cr}}(i, y_{1850-2019})}{A_{\text{HYDE-cr}}(i, y_{2000})} \qquad (16)$$

$$R_{\text{HYDE-past}}(i, y_{1850-2019}) = \frac{A_{\text{HYDE-past}}(i, y_{1850-2019})}{A_{\text{HYDE-past}}(i, y_{2000})} \qquad (17)$$

where $y_{2000}$ is the year 2000; $y_{1850-2019}$ is the year (in the period 1850–2019).

We then applied these normalised (temporal) values of cropland and pasture area of Eqs. (16–17) to the respective land use area of Ramankutty et al.³⁴ in 2000, as given in Eqs. (18) and (19). By doing so, we maintained the spatial distribution of land use areas as of Ramankutty et al.³⁴, while considering the annual temporal developments from HYDE. Finally, we derived agriculture area ($A_{\text{agri}}$ (ha)) by combining our reconstructed gridded cropland and pasture area ($A_{\text{cr}}$ (ha) and $A_{\text{past}}$ (ha), respectively) for the time period 1850–2019, as given in Eq. (20).

$$A_{\text{cr}}(i, y_{1850-2019}) = A_{\text{Ramankutty-crop}}(i, y_{2000}) \times R_{\text{HYDE-cr}}(i, y_{1850-2019}) \qquad (18)$$

$$A_{\text{past}}(i, y_{1850-2019}) = A_{\text{Ramankutty-past}}(i, y_{2000}) \times R_{\text{HYDE-past}}(i, y_{1850-2019}) \qquad (19)$$

$$A_{\text{agri}}(i, y_{1850-2019}) = A_{\text{cr}}(i, y_{1850-2019}) + A_{\text{past}}(i, y_{1850-2019}) \qquad (20)$$

Harmonization of cropland, pasture and agriculture area at country-level with FAOSTAT. We used the FAOSTAT dataset as reference database for country-level information, due to its consistent availability for the period 1961–2019 and global coverage across a range of variables required to estimate the N surplus. To match our estimate of cropland and pasture area, we first derived “Cropland” and “Land and permanent meadows and pastures” from FAOSTAT for the time period 1961–2019. Then, for each year and country, we calculated country-level differences (as ratios) for cropland and pasture ($R_{A_{\text{cr}}}$ (-) and $R_{A_{\text{past}}}$ (-), respectively) between those given in FAOSTAT ($A_{\text{FAO-cr}}$ (ha) and $A_{\text{FAO-past}}$ (ha), respectively) and ones estimated in our study ($A_{\text{cr}}$ (ha) and $A_{\text{past}}$ (ha)), as summarized in Eqs. (21–22):

$$R_{A_{\text{cr}}}(u, y_{1961-2019}) = \frac{A_{\text{FAO-cr}}(u, y_{1961-2019})}{\sum_{i=1}^{n_u} A_{\text{cr}}(i, y_{1961-2019})} \qquad (21)$$

$$R_{A_{\text{past}}}(u, y_{1961-2019}) = \frac{A_{\text{FAO-past}}(u, y_{1961-2019})}{\sum_{i=1}^{n_u} A_{\text{past}}(i, y_{1961-2019})} \qquad (22)$$

where $y_{1961}$ is the year 1961; $y_{1961-2019}$ is the year (in the period 1961–2019). u is country; $n_u$ is the number of grid cell in the u-th country.

Then, we applied these calculated ratio to our gridded estimates of cropland and pasture areas to ensure that the FAOSTAT country totals are maintained for the period 1961–2019. The resulting gridded cropland and pasture areas ($C_{A_{\text{cr}}}$ (ha) and $C_{A_{\text{past}}}$ (ha), respectively) can be given in Eqs. (23–24) as:

$$C_{A_{\text{cr}}}(i, y_{1961-2019}) = R_{A_{\text{cr}}}(u, y_{1961-2019}) \times A_{\text{cr}}(i, y_{1961-2019}) \qquad (23)$$

$$C_{A_{\text{past}}}(i, y_{1961-2019}) = R_{A_{\text{past}}}(u, y_{1961-2019}) \times A_{\text{past}}(i, y_{1961-2019}) \qquad (24)$$

As FAOSTAT does not provide estimates before 1961, we applied the same ratio as of 1961 for the time period 1850–1960 as given in Eqs. (25–26):

$$C_{A_{\text{cr}}}(i, y_{1850-1960}) = R_{A_{\text{cr}}}(u, y_{1961}) \times A_{\text{cr}}(i, y_{1850-1960}) \qquad (25)$$

$$C_{A_{\text{past}}}(i, y_{1850-1960}) = R_{A_{\text{past}}}(u, y_{1961}) \times A_{\text{past}}(i, y_{1850-1960}) \qquad (26)$$

For the countries for which FAOSTAT data are missing before 1992, such as Croatia, Estonia, Latvia, Lithuania and Slovenia, we applied the same ratio as of the year 1992 for the period 1850–1991. Furthermore, for the countries like Belgium and Luxembourg, Czech Republic and Slovakia, the FAOSTAT maintains single (combined) values of reported variables for the past records (prior to 1993 for Czechoslovakia and before 2000 for Belgium-Luxembourg). In our estimation of country-specific ratios we took care of these details, and accordingly applied a single ratio factor for the adjoining countries and records. Finally, the agricultural area was derived as a sum of the FAOTSTAT harmonized cropland and pasture area as mentioned in Eq. (27):
$$C_{A_{\text{agri}}}(i, y_{1850-2019}) = C_{A_{cr}}(i, y_{1850-2019}) + C_{A_{\text{past}}}(i, y_{1850-2019}) \quad (27)$$

where $y_{1961}$ is the year 1961; $y_{1850-1960}$ is the year (in the period 1850–1960).

Finally, we checked for a physical consistency of grid-level reconstructed area estimates i.e., the exceedance of the gridded agriculture area of Eq. (27) ($C_{A_{\text{agri}}}(ha)$) beyond the physical area of a given grid cell. This condition, although appeared seldom, was unavoidable due to differences in origin of different databases (e.g., FAOSTAT, HYDE, Ramankutty et al.³⁴) used in our reconstruction. To maintain the physical consistency, we equally redistributed the agricultural area in excess to other grid cells in a way that the agricultural area in a grid cell is not greater than the area of a grid cell. In a similar way, we corrected the gridded cropland areas ($C_{A_{cr}}(ha)$) and pasture areas ($C_{A_{\text{past}}}(ha)$) to guarantee that FAOSTAT country totals are satisfied and that the sum of cropland area and pasture area does not exceed the total agricultural area in a given grid cell.

**Reconstruction of non-agriculture area.** We derived the non-agricultural area in a given grid cell as the area $A_{\text{NonAgri}}(ha)$ that is left after allocation of cropland and pasture area, as given in Eq. (28):

$$A_{\text{NonAgri}}(i, y_{1850-2019}) = 1 - [C_{A_{cr}}(i, y_{1850-2019}) + C_{A_{\text{past}}}(i, y_{1850-2019})] \quad (28)$$

Afterwards, we utilized the classification of land cover categories from global land cover (GLC)³⁹ that is available at a spatial resolution of 300 m and for the year 2000. GLC includes 23 land cover classes that we grouped into 5 categories namely, cropland, semi-natural-vegetation (tree, shrub-land, herbaceous cover), forest (broad-leaved, evergreen and deciduous forest), non-vegetation (bare areas, water bodies) and urban area. We calculated the fractions of forest, semi-natural-vegetation, non-vegetation and urban area (ForGLC NatVegGLC NonVegGLC UrbanGLC, -) respectively) from GLC land cover within our target grid cell of 5 arcmin, and their sum which corresponds to the total non-agricultural fraction (NonAgriGLC, -). We then estimated the relative proportion of forest, semi-natural vegetation, non-vegetation and urban area $fr_{ForGLC}$, $fr_{NatVegGLC}$, $fr_{NonVegGLC}$, $fr_{UrbanGLC}$, - respectively) within the total non-agricultural fraction as given from Eqs. (29) to (32):

$$fr_{ForGLC}(i, y_{2000}) = \frac{ForGLC(i, y_{2000})}{NonAgriGLC(i, y_{2000})} \quad (29)$$

$$fr_{NatVegGLC}(i, y_{2000}) = \frac{NatVegGLC(i, y_{2000})}{NonAgriGLC(i, y_{2000})} \quad (30)$$

$$fr_{UrbanGLC}(i, y_{2000}) = \frac{UrbanGLC(i, y_{2000})}{NonAgriGLC(i, y_{2000})} \quad (31)$$

$$fr_{NonVegGLC}(i, y_{2000}) = \frac{NonVegGLC(i, y_{2000})}{NonAgriGLC(i, y_{2000})} \quad (32)$$

To infer the annual development of forest, semi-natural vegetation, non-vegetation and urban land cover areas ($A_{For}$, $A_{NatVeg}$, $A_{NonVeg}$ and $A_{Urban}$, $ha$), respectively), we multiplied the corresponding relative fractions (Eqs. (29) to (32)) with the gridded non-agricultural area ($A_{NonAgri}$, $ha$) of Eq. (28), as summarised in Eqs. (33) to (36):

$$A_{For}(i, y_{1850-2019}) = fr_{For}(i, y_{2000}) \times A_{NonAgri}(i, y_{1850-2019}) \quad (33)$$

$$A_{NatVeg}(i, y_{1850-2019}) = fr_{NatVeg}(i, y_{2000}) \times A_{NonAgri}(i, y_{1850-2019}) \quad (34)$$

$$A_{NonVeg}(i, y_{1850-2019}) = fr_{NonVeg}(i, y_{2000}) \times A_{NonAgri}(i, y_{1850-2019}) \quad (35)$$

$$A_{Urban}(i, y_{1850-2019}) = fr_{Urban}(i, y_{2000}) \times A_{NonAgri}(i, y_{1850-2019}) \quad (36)$$

**Crop specific harvested area.** Non-fodder crops. We obtained gridded crop specific harvested areas from Monfreda et al.³⁵ which is available for 175 different crop types and that correspond to the representative year of 2000. Monfreda et al.³⁵ generated crop harvested area by combining information from surveys and agricultural census for around 206 countries. In addition to national statistics, a substantial information in these databases is derived from Agro-MAPS which is a project of combined efforts from FAOSTAT, the International Food Policy Research Institute (IFPRI), and the Center for Sustainability and the Global Environment (SAGE)³⁵. Among these 175 crops, we selected 17 major crops for which fertilizer application rates are provided⁴⁰ and that have either a large production across Europe or a N content greater than 1 kg of N tonne⁻¹ of product. The selected crops are listed in Table 2.
<table><caption>**Table 2.** N content and biological N fixation rates for crops and other land types included in this study. The values for N content and BNF rate were derived based on estimates reported in previous studies<sup>8,18,24,68</sup>.</caption><thead><tr><th rowspan="2">Crop/Land type</th><th>N content</th><th>BNF rate</th></tr><tr><th>kg of N t<sup>-1</sup> of product</th><th>kg of N ha<sup>-1</sup> of area</th></tr></thead><tbody><tr><td colspan="3"><strong>Non-fodder crops</strong></td></tr><tr><td>Wheat</td><td>20<sup>8,18</sup></td><td>—</td></tr><tr><td>Maize</td><td>13<sup>8,18</sup></td><td>—</td></tr><tr><td>Rice</td><td>12<sup>8,24</sup></td><td>—</td></tr><tr><td>Soybean</td><td>45<sup>8,18</sup></td><td>80<sup>18,68</sup></td></tr><tr><td>Sunflower seed</td><td>23<sup>18,24</sup></td><td>—</td></tr><tr><td>Millet</td><td>24<sup>18,24</sup></td><td>—</td></tr><tr><td>Sorghum</td><td>14<sup>24</sup></td><td>—</td></tr><tr><td>Rye</td><td>22<sup>18,24</sup></td><td>—</td></tr><tr><td>Barley</td><td>18<sup>18,24</sup></td><td>—</td></tr><tr><td>Oats</td><td>20<sup>18,24</sup></td><td>—</td></tr><tr><td>Triticale, Buckwheat</td><td>16<sup>24</sup></td><td>—</td></tr><tr><td>Rapeseed, Sesame</td><td>35<sup>24</sup></td><td>—</td></tr><tr><td>Pulses</td><td>40<sup>8,18</sup></td><td>50<sup>8,18</sup></td></tr><tr><td>Sugar-beet</td><td>1.5<sup>18,24</sup></td><td>—</td></tr><tr><td>Potato</td><td>3<sup>18,24</sup></td><td>—</td></tr><tr><td colspan="3"><strong>Fodder crops</strong></td></tr><tr><td>Temporary grassland, lucerne, green maize, root crops etc.</td><td>—<sup>21</sup></td><td>—<sup>21</sup></td></tr><tr><td colspan="3"><strong>Forest</strong></td></tr><tr><td>Temperate</td><td>—</td><td>16<sup>54</sup></td></tr><tr><td>Boreal</td><td>—</td><td>1.77<sup>54</sup></td></tr><tr><td><strong>Pasture</strong></td><td>—</td><td>5<sup>8</sup></td></tr><tr><td><strong>Semi-natural vegetation</strong></td><td>—</td><td>2.7<sup>54</sup></td></tr></tbody></table>

To construct annual estimates of crop specific harvested area ($A_{\text{crops}}$ (ha)), we first applied temporal dynamics of cropland area ($C_{A_{\text{cr}}}$ (ha)) derived in Eqs. (23) and (25) with respect to the gridded dataset from Monfreda et al.<sup>35</sup> ($A_{\text{crops}_{\text{Monfreda}}}$ (ha)); as given in Eq. (37). Consequently, the temporal dynamics in our estimates are referenced to the year 2000 to maintain the consistency with the Monfreda et al.<sup>35</sup> database and the overall methodology applied in the reconstruction of cropland/agricultural areas (as described above).

$$A_{\text{crops}}(i, c, y_{1850-2019}) = A_{\text{crops}_{\text{Monfreda}}}(i, c, y_{2000}) \times \frac{C_{A_{\text{cr}}}(i, y_{1850-2019})}{C_{A_{\text{cr}}}(i, y_{2000})} \quad (37)$$

where c is crop, $y_{1850-2019}$ is the year (in the period 1850–2019).

We then ensured that the calculated crop specific harvested areas match estimates from FAOSTAT at a country level. To do so, we derived the country level ratio ($R_A$ (-)) between crop harvested area from FAOSTAT ($A_{\text{crops}_{\text{FAO}}}$ (ha)) and crop area provided in this study, as given in Eq. (38):

$$R_A(u, c, y_{1961-2019}) = \frac{A_{\text{crops}_{\text{FAO}}}(u, c, y_{1961-2019})}{\sum_{i=1}^{n_u} A_{\text{crops}}(u, c, y_{1961-2019})} \quad (38)$$

The calculated ratio was then applied to the estimates of gridded annual crop specific harvested area ($A_{\text{crops}}$ (ha)) to obtain the crop harvested area harmonized at country level with FAOSTAT ($C_{A_{\text{crops}}}$ (ha)), as mentioned in Eq. (39). Since FAOSTAT data only starts in 1961, all years prior to 1961 have been rescaled by the scaling factor (ratio) of 1961 as given in Eq. (40):

$$C_{A_{\text{crops}}}(i, c, y_{1961-2019}) = A_{\text{crops}}(i, c, y_{1961-2019}) \times R_A(u, c, y_{1961-2019}) \quad (39)$$

$$C_{A_{\text{crops}}}(i, c, y_{1850-1960}) = A_{\text{crops}}(i, c, y_{1850-1960}) \times R_A(u, c, y_{1960}) \quad (40)$$

**Fodder crops.** Since FAOSTAT does not provide data on fodder crops, we used for this the crop areas from Einarsson et al.<sup>21</sup>. This dataset is available for 26 European countries during 1961–2019. The dataset is given for 6 fodder crop categories (Table 2) namely: temporary grassland, lucerne, other leguminous plants, green maize, plants harvested from arable land and other root crops (forage beet, turnip etc.). For the temporal development of fodder crop area during 1850–1960, we utilized the temporal variability of the reconstructed cropland area
(see Eqs. (23) and (25)). To do so, we first derived country-level ratio ($R_{A_{fodder}}$ (-)) of fodder crop area ($A_{fodder}$ ($ha$)) to the cropland area ($C_{A_{cr}}$ ($ha$)) in the year 1961 as given in Eq. (41). Later, we applied the country-level ratio of the year 1960, to the cropland year during 1850–1960, as given in Eq. (42).

$$R_{A_{fodder}}(u, c, y_{1961}) = \frac{A_{fodder}(u, c, y_{1961})}{\sum_{i=1}^{n_u} C_{A_{cr}}(i, y_{1961})} \quad (41)$$

$$A_{fodder}(u, c, y_{1850-1960}) = R_{A_{fodder}}(u, c, y_{1961}) \times \sum_{i=1}^{n_u} C_{A_{cr}}(i, y_{1850-1960}) \quad (42)$$

where c is fodder crop, $y_{1961}$ is the year 1961, $y_{1850-1960}$ is the year (in the period 1850–1960).

Einarsson et al.²¹ does not provide fodder crop area before 1992 for some countries such as Croatia, Estonia, Latvia, Lithuania, Slovenia, Czechia, Slovakia, Belgium and Luxembourg. To fill the data gaps in fodder crop area, following consideration were made:

* Where possible, we extrapolated country level ratio ($R_{A_{fodder}}$ (-)) from neighbouring countries with similar climate and geographical conditions. For instance, the country level ratio ($R_{A_{fodder}}$ (-)) of Austria, Sweden, Romania, Greece was used for Switzerland, Norway, Ukraine, Turkey, respectively.

* For the other missing countries that are located in Eastern Europe (e.g. Croatia, Estonia), we used an aggregated country-level ratio ($R_{A_{fodder}}$ (-)) from other available countries in Eastern Europe namely Bulgaria, Romania, Hungary for the year 1992. This ratio was then applied to the cropland area of the missing countries during 1850–1991.

We thus completed the time series of fodder crop area for the period 1850–2019 at country-level, that we downscaled to gridded level using the fodder specific crop area from Monfreda et al.³⁵, similar to the approach described above for the reconstruction of non-fodder crop areas.

**N inputs.** N inputs in our estimates consisted of mineral fertilizer, animal manure, atmospheric deposition and biological fixation to agriculture and non-agriculture soil at gridded level from 1850 to 2019.

**Fertilizer.** The amount of fertilizers applied to croplands and pastures is generally assessed using application rates for the different crop types and pastures¹,³⁰. In particular, a recent study suggests that a large share of fertilizers are applied to permanent pastures in many European countries²¹. However, differences in the fractions of fertilizer applied to croplands and pastures have been reported as being a major source of uncertainty⁵. To account for these uncertainties, we derived two gridded estimates of fertilizer application based on two different sets of application rates for croplands and pastures. These two sets of rates were then used to downscale the country-level fertilizer amounts applied to soil at a gridded level.

**Country-level fertilizer applied to soil during 1850–2019.** During 1961 to 2019, we derived country-level data of fertilizer applied to agricultural soils from FAOSTAT⁴¹. For a few countries, namely Croatia, Estonia, Latvia, Lithuania, Belarus, Ukraine, no data are provided before 1992. For these countries, we estimated the N fertilizer application ($N \text{fer}_{Missing_{EastEU}}$ ($kg \text{yr}^{-1}$)) during the period 1961–1991 by applying the temporal dynamics of Eastern European countries that have reported data i.e., Bulgaria, Hungary and Czechoslovakia, as given in Eq. (43):

$$N \text{fer}_{Missing_{EastEU}}(u, y_{1961-1991}) = \frac{N \text{fer}_{EastEU}(y_{1961-1991})}{N \text{fer}_{EastEU}(y_{1992})} \times N \text{fer}_{Missing_{EastEU}}(u, y_{1992}) \quad (43)$$

where $N \text{fer}_{EastEU}$ ($kg \text{yr}^{-1}$) is the total fertilizer application calculated over Bulgaria, Hungary and Czechoslovakia. Moreover, Belgium and Luxembourg are reported in FAOSTAT as a single country during 1961–1999, while separate country estimates are provided between 2000–2019. Similarly, Czechia and Slovakia is reported in a single database of Czechoslovakia between 1961–1992. Before 2000 for Belgium and Luxembourg and before 1993 for Czechia and Slovakia, we applied the dynamics of Belgium-Luxembourg and former Czechoslovakia, respectively.

During 1850–1960, we lacked country-level information from FAOSTAT. During 1850–1920, we assumed that there is no fertilizer applied to agricultural soil based on evidence from previous studies¹,⁴². Regarding the time period of 1920–1960, we utilized the temporal dynamics from Holland et al.⁴² that provides global estimates of fertilizer production during 1925–1960 based on different data sources⁴³,⁴⁴. We filled the values of global fertilizer production for the years between 1920 and 1924 using linear interpolation with no (or zero) fertilizer in the year 1920 and Holland et al.⁴² estimates for the year 1925. The global temporal dynamic was then applied to all countries in the study domain between 1920–1960 with respect to their corresponding estimates for the year 1961. We termed the completed annual country-level fertilizer amount as $N \text{fer}_{soil}(u, y_{1850-2019})$.

**Distribution of fertilizer to croplands and pastures.** To downscale the annual country-level fertilizer amount to croplands and pastures at grid level, we utilized two different sets of application rates to account for the
uncertainties in the gridded spatial pattern of fertilizer application within a given country. We first derived country-level fertilizer application rates for individual crops and grassland based on information given by the International Fertilizer Industry Association (IFA; https://www.ifastat.org), which are then adjusted using two different approaches. We provide further details in the paragraphs below.

We obtained country-level fertilizer use per crop types and grassland ($N fer_{crops}$ and $N fer_{grass}$, respectively ($kg yr^{-1}$)) from IFA in 2014–2015$^{40}$. IFA provides country-level total amount of N fertilizer use for 13 crop groups, i.e., wheat, rice, maize (both grain and silage), other cereals (barley, sorghum, oats, rye, triticale, millet, etc.), soybean, oil palm, other oil seeds (rapeseed, mustard, sunflower, groundnut, coconut, etc), fiber crops (cotton, flax, hemp, jute, etc), sugar crops (sugar cane and sugar beet), roots and tubers, fruits, vegetables, grassland (pastures for hay, silage and grazing, and temporary grasslands) and residual (pulses, nut trees, rubber, cocoa, coffee, tea, tobacco, etc.). In Europe, IFA data does not contain the N fertilizer use for each individual country but it reports the total value for EU-28, as well as values for Ukraine, Russia and Belarus. We combined fertilizer use amount from IFA and harvested areas of crops and grassland from FAOSTAT$^{38}$ ($A_{crops}$ and $A_{grass}$ ($ha$), respectively) to estimate the respective application rates. The area of grassland includes temporary grasslands (“land under temporary meadows and pastures” in FAOSTAT) and pastures (“land under permanent meadows and pastures” in FAOSTAT). We thus derived country-level fertilizer application rates for individual non-fodder crops ($N fer_{cropsRate}$ ($kg ha^{-1}$)) of crop harvested areas ($yr^{-1}$) and grassland ($N fer_{grassRate}$ ($kg ha^{-1}$)) of grassland areas ($yr^{-1}$), as summarized in Eqs. (44–45):

$$N fer_{cropsRate}(u, c, y_{2015}) = \frac{N fer_{crops}(u, c, y_{2015})}{A_{crops}(u, c, y_{2015})} \quad (44)$$

$$N fer_{grassRate}(u, y_{2015}) = \frac{N fer_{grass}(u, y_{2015})}{A_{grass}(u, y_{2015})} \quad (45)$$

For fodder crops and pastures, we adopted the fertilizer application rates for the grassland (Eq. (45)). For the rest of the countries in our domain that are not the part of the IFA dataset, we assumed the same fertilizer application rates for individual crops and pastures as of EU-28.

Afterwards, to inform the spatial variability, we multiplied the country-level fertilizer application rates for non-fodder crops ($N fer_{cropsRate}$) and grassland ($N fer_{grassRate}$) by the gridded areas of non-fodder crops, pastures and fodder crops ($C_{A_{crops}}$, $C_{A_{past}}$ and $A_{fodder}$ respectively, ($ha$)) during 1850–2019. Consequently, we obtained annual fertilizer application amount for each crop type (non-fodder and fodder), pastures, and total ($N fer_{crops}$, $N fer_{fodder}$, $N fer_{past}$, $N fer_{soil}$, respectively ($kg yr^{-1}$)) for each grid cell, as summarised in Eqs. (46–49):

$$N fer_{crops}(i, c, y_{1850-2019}) = N fer_{cropsRate}(u, c, y_{2015}) \times C_{A_{crops}}(i, c, y_{1850-2019}) \quad (46)$$

$$N fer_{past}(i, y_{1850-2019}) = N fer_{grassRate}(u, y_{2015}) \times C_{A_{past}}(i, y_{1850-2019}) \quad (47)$$

$$N fer_{fodder}(i, c, y_{1850-2019}) = N fer_{grassRate}(u, y_{2015}) \times A_{fodder}(i, c, y_{1850-2019}) \quad (48)$$

$$N fer_{soil}(i, y_{1850-2019}) = \sum_{i=1}^{n_c} N fer_{crops}(i, c, y_{1850-2019}) + N fer_{past}(i, y_{1850-2019}) + \sum_{i=1}^{n_c} N fer_{fodder}(i, c, y_{1850-2019}) \quad (49)$$

In the next step, we adjusted the fertilizer application amount (calculated in Eq. (49)) to ensure that it is consistent with the country-level fertilizer applied to soil during 1850–2019 reconstructed in previous steps ($N fer_{soil}$ ($kg yr^{-1}$)). For this purpose, we first derived an adjustment factor (as ratio) of the country-level fertilizer amount to gridded fertilizer amount aggregated to country level. The calculated ratio was then multiplied by the fertilizer application rates of individual crops and grassland ($N fer_{cropsRate}$ ($kg ha^{-1}$) of crop harvested areas ($yr^{-1}$)) and $N fer_{grassRate}$ ($kg ha^{-1}$) of grassland areas ($yr^{-1}$), respectively) to derive adjusted fertilizer application rates for crops and grassland ($C_{N fer_{cropsRate}}$ ($kg ha^{-1}$) of crop harvested areas ($yr^{-1}$)) and $C_{N fer_{grassRate}}$ ($kg ha^{-1}$) of grassland areas ($yr^{-1}$), respectively), as given from Eqs. (50–51):

$$C_{N fer_{cropsRate}}(u, c, y_{1850-2019}) = N fer_{cropsRate}(u, c, y_{2015}) \times \frac{N fer_{soil}(u, y_{1850-2019})}{\sum_{i=1}^{n_u} N fer_{soil}(i, y_{1850-2019})} \quad (50)$$

$$C_{N fer_{grassRate}}(u, y_{1850-2019}) = N fer_{grassRate}(u, y_{2015}) \times \frac{N fer_{soil}(u, y_{1850-2019})}{\sum_{i=1}^{n_u} N fer_{soil}(i, y_{1850-2019})} \quad (51)$$
Now to determine the partitioning of fertilizer application amount to croplands and pastures, we considered two different sets of application rates for croplands and pastures to account for (methodological) uncertainties. In the first approach, we adopted the application rates based on IFA rates along with their subsequent adjustments as detailed above in Eqs. (50) and (51). In the second approach, we further adjusted these rates so that the proportions of fertilizer applied to croplands and pastures are consistent with the (partitioning) information provided by Einarsson et al.²¹ for the period 1961–2019. For the countries with incomplete record in Einarsson et al.²¹, we used the proportions from either neighbouring countries or the average value of countries in Eastern Europe to fill the missing years, similar to the procedure we adopted in section “Crop specific harvested area-Fodder crops”.

Finally, from the two sets of fertilizer application rates, we estimated the gridded fertilizer amount applied to croplands and pastures (N fer<sub>cr</sub> (kg yr<sup>-1</sup>) and N fer<sub>past</sub> (kg yr<sup>-1</sup>), respectively) using the gridded areas of non-fodder crops, fodder crops and pastures. For instance, for the first approach (using the application rates of Eqs. (50) and (51)), the equations can be summarized as:

$$ N fer_{cr}(i, y_{1850-2019}) = \sum_{i=1}^{n_c} (C_{N fer_{cropsRate}}(u, c, y_{1850-2019}) \times C_{A_{crops}}(i, c, y_{1850-2019})) + \sum_{i=1}^{n_c} (C_{N fer_{grassRate}}(u, y_{1850-2019}) \times A_{fodder}(i, c, y_{1850-2019})) \quad (52) $$

$$ N fer_{past}(i, y_{1850-2019}) = C_{N fer_{grassRate}}(u, y_{1850-2019}) \times C_{A_{past}}(i, y_{1850-2019}) \quad (53) $$

For the second approach, we replaced the fertilizer application rates that are used in Eqs. (52) and (53) (C<sub>N fer<sub>cropsRate</sub></sub> and C<sub>N fer<sub>grassRate</sub></sub>) by the adjusted rates based on information provided by Einarsson et al.²¹.

Within each approach, the total gridded fertilizer amount applied to soil was then derived by summing the fertilizer applied to croplands and pastures. We thus obtained two different datasets of fertilizer application that reflect the methodological uncertainties in the gridded spatial pattern within a given country, while the country-level values are the same in the both datasets.

**Animal manure.** The amount of N excreted by the livestock (also called manure production) is typically calculated using N excretion rates and livestock counts. In addition, it can have different fates i.e., it can be left on pasture or it can be collected, stored and then applied to soils (cropland and pasture). The differences in livestock N excretion amounts obtained using different N excretion rates and the distribution of manure production between cropland and pasture has been recognized as another major source of uncertainty in N budget estimations⁵,²⁷,²⁸.

To account for these uncertainties, we utilized information from two datasets of animal manure (i.e., FAOSTAT⁴⁵ and Einarsson et al.²¹) available at country-level during 1961–2019. Both datasets differ in estimating the total N excretion from livestock category due to different data input sources and methodology. For example, FAOSTAT calculates the total amount of N excreted by the livestock categories by multiplying the country-level livestock counts by regional-level values of typical animal mass and N excretion rates from the Intergovernmental Panel on Climate Change (IPCC)⁴⁶. Einarsson et al.²¹, on the other hand, follows the approach of Lassaletta et al.⁴⁷ in which the N excretion rates within a given livestock category are assumed to be proportional to their slaughter weight. Hence the N excretion rates are derived by multiplying the slaughtered animal weights by the ratio of N excretion rate to slaughtered animal weight for each country and year. The derived N excretion rates are then multiplied by the livestock counts to calculate total N excretion in each country and year. We included both these databases in our analysis. Note that for the missing estimates (for countries and time-periods) in the database from Einarsson et al.²¹, we adapted the same methodology as explained in section “Crop specific harvested area-Fodder crops”. We derived country-level data of manure applied to soils (Man<sub>app</sub> (kg yr<sup>-1</sup>)) and manure left on pastures (Man<sub>left past</sub> (kg yr<sup>-1</sup>)) for the period 1961–2019 from FAOSTAT⁴⁸ and Einarsson et al.²¹. As we focus on the soil N-budget, the manure applied to soils in both databases excludes losses of reduced forms of N via volatilization occurring during manure management and storage. Subsequent deposition of these reduced forms of N is accounted for in the budget in the estimates of atmospheric N deposition.

To downscale the country-level values to gridded level, we utilized gridded information of manure production from Zhang et al.²⁹ that was derived using the spatial distribution of livestock counts from the Global Livestock Impact Mapping System (GLIMS)⁴⁹ and N excretion coefficients from the IPCC at a 5 arcmin spatial resolution for the time period 1860–2014. Zhang et al.²⁹ does not give estimate of gridded manure production before 1860 and after 2015. We assumed the same amount of manure production for the time period 1850–1859 as of the year 1860 and for the time period 2015–2019 as of the year 2014. We combined the gridded manure production (Man<sub>prod</sub> (kg ha<sup>-1</sup>) of grid area yr<sup>-1</sup>)) of Zhang et al.²⁹ and the grid cell area (A<sub>grid</sub> (ha)) to calculate the country-level ratio of manure applied to soil to manure production (R<sub>app,grid/prod</sub> (-)) and manure left on pastures to manure production (R<sub>left past/prod</sub> (-)), given below in the Eqs. (54–55). Note that these ratios were estimated for each of two manure databases²¹,⁴⁵, separately.
$$R_{app_{soil}/prod}(u, y_{1961-2019}) = \frac{Man_{app_{soil}}(u, y_{1961-2019})}{\sum_{i=1}^{n_u} Man_{prod}(i, y_{1961-2019}) \times A_{grid}(i)} \quad (54)$$

$$R_{left_{past}/prod}(u, y_{1961-2019}) = \frac{Man_{left_{past}}(u, y_{1961-2019})}{\sum_{i=1}^{n_u} Man_{prod}(i, y_{1961-2019}) \times A_{grid}(i)} \quad (55)$$

Next, we estimated the gridded rates of manure applied to soil and left on pastures (Man$_{app_{soil}$}$ and Man$_{left_{past}$}$, respectively, both are given in kg ha$^{-1}$ of grid area yr$^{-1}$) by multiplying the calculated ratios from Eqs. (54) and (55) by the gridded manure production from Zhang et al.²⁹ (Man$_{prod}$ (kg ha$^{-1}$) of grid area yr$^{-1}$)), as presented in Eqs. (56–57).

$$Man_{app_{soil}}(i, y_{1961-2019}) = R_{app_{soil}/prod}(u, y_{1961-2019}) \times Man_{prod}(i, y_{1961-2019}) \quad (56)$$

$$Man_{left_{past}}(i, y_{1961-2019}) = R_{left_{past}/prod}(u, y_{1961-2019}) \times Man_{prod}(i, y_{1961-2019}) \quad (57)$$

For the time period during 1850–1960, the same ratio of the year 1961 was applied to the manure production database, as given in Eqs. (58) and (59).

$$Man_{app_{soil}}(i, y_{1850-1960}) = R_{app_{soil}/prod}(u, y_{1961}) \times Man_{prod}(i, y_{1850-1960}) \quad (58)$$

$$Man_{left_{past}}(i, y_{1850-1960}) = R_{left_{past}/prod}(u, y_{1961}) \times Man_{prod}(i, y_{1850-1960}) \quad (59)$$

Finally, we checked the consistency of gridded manure application rates to fall within a maximum permissible limit following Bouwman et al.⁵⁰ for European countries (170–250 kg ha$^{-1}$ of agricultural area yr$^{-1}$). In grid cells and years, where this maximum limit was crossed, the extra manure was redistributed to the neighbouring agricultural cells.

*Distribution of manure applied between cropland and pasture.* To distribute manure applied to soil (calculated in Eqs. (56) and (58) from both FAOSTAT and Einarsson et al.²¹ data) between cropland and pasture, we used two different methodologies to account for uncertainties. First, following a similar approach as Tian et al.⁵¹, we assumed the uniform manure application rates between cropland and pasture within a given grid cell. Thus, manure applied to cropland (Man$_{cr}$ (kg yr$^{-1}$)) was manure applied to soil (Man$_{app/soil}$ (kg ha$^{-1}$ yr$^{-1}$)) times crop-land area (C$_{Acr}$ (ha)) in a given grid cell according to Eq. (60). Similarly, manure applied to pastures (Man$_{app_{past}}$ (kg yr$^{-1}$)) was manure applied to soil (Man$_{app/soil}$ (kg ha$^{-1}$ yr$^{-1}$)) times pasture area (C$_{A_{past}}$ (ha)) in a given grid cell, as given in Eq. (61). The final manure amount to pastures (Man$_{past}$ (kg yr$^{-1}$)) was derived as a sum of manure applied and left on pasture as expressed in Eq. (62).

$$Man_{cr}(i, y_{1850-2019}) = Man_{app/soil}(i, y_{1850-2019}) \times C_{A_{cr}}(i, y_{1850-2019}) \quad (60)$$

$$Man_{app_{past}}(i, y_{1850-2019}) = Man_{app/soil}(i, y_{1850-2019}) \times C_{A_{past}}(i, y_{1850-2019}) \quad (61)$$

$$Man_{past}(i, y_{1850-2019}) = Man_{app_{past}}(i, y_{1850-2019}) + Man_{left_{past}}(i, y_{1850-2019}) \quad (62)$$

In the second manure application approach, we distributed the manure applied to soil (calculated in Eqs. (56) and (58) from both FAOSTAT and Einarsson et al.²¹ data) according to the country-level proportion of manure applied to cropland and pasture derived from Einarsson et al.²¹. The latter study²¹ allocated stored manure to cropland and pasture using country-specific expert estimates of the share of manure applied to cropland and pasture based on national statistics collected for different animal categories for given countries. We derived country-level, annual ratio of manure applied to cropland (and manure applied to pasture) to total manure applied to soil. We then adjusted the gridded amounts of manure applied to cropland and pasture (Eqs. (60) and (62), respectively) to comply with Einarsson et al.²¹ estimated ratios.

Overall, by utilizing two different input data sources of manure applied to soil and left on pasture (FAOSTAT and Einarsson et al.²¹) and two different methodologies to distribute manure applied to soil between croplands and pastures, we reconstructed four different estimates of gridded manure in our database.

*Atmospheric deposition.* Both reduced and oxidized forms of N from the atmosphere fall back to terrestrial and aquatic ecosystem through either wet deposition (i.e, by precipitation or snow) or dry deposition by process termed N deposition⁵². Atmospheric N deposition in this study was quantified using the National Center for Atmospheric Research (NCAR), Chemistry–Climate Model Initiative (CCMI) N-deposition dataset, which is part of the input datasets for Model Inter comparison Projects (input4MIPS)⁵³. The dataset is given at monthly time step from 1850 to 2014 at a spatial resolution of 1.9° latitude by 2.5° longitude. We downscaled the data to
the required grid cell size of 5 arcmin resolution using nearest neighbor interpolation and also aggregated to annual time step. To get N deposition for cropland, pasture, forest, semi-natural vegetation, non-vegetation and urban land use, we multiplied the total N deposition with the proportion of respective land area per grid cell.

**Biological nitrogen fixation (BNF).** BNF is a natural process of transforming atmospheric non-reactive nitrogen ($N_2$) to its reactive form (ammonia ($NH_3$) or nitrates) through microbial activity⁵⁴,⁵⁵. In this study, N fixation over cropland, pasture, forest and semi-natural vegetation ($BNF_{cr}$, $BNF_{pasts}$, $BNF_{For}$ and $BNF_{NatVegs}$ respectively) was derived by multiplying the respective areas ($C_{A_{crops}}$, $C_{A_{past}}$, $A_{For}$ and $A_{NatVeg}$ ($ha$), respectively) with their N fixation rates ($BNF_{crops-Rate}$, $BNF_{past-Rate}$, $BNF_{For-Rate}$ and $BNF_{NatVeg-Rate}$ ($kg \ ha^{-1} \ yr^{-1}$), respectively); as summarised in Eqs. (63–66).

$$BNF_{cr}(i, y_{1850-2019}) = \sum_{i=1}^{n_c} (C_{A_{crops}}(i, c, y_{1850-2019}) \times BNF_{crops-Rate}(c)) \quad (63)$$

$$BNF_{past}(i, y_{1850-2019}) = C_{A_{past}}(i, y_{1850-2019}) \times BNF_{past-Rate} \quad (64)$$

$$BNF_{For}(i, y_{1850-2019}) = A_{For}(i, y_{1850-2019}) \times BNF_{For-Rate} \quad (65)$$

$$BNF_{NatVeg}(i, y_{1850-2019}) = A_{NatVeg}(i, y_{1850-2019}) \times BNF_{NatVeg-Rate} \quad (66)$$

N fixation rates are derived for different crops and land use categories from previous studies⁸,²¹,⁵⁴. For fodder crops, we calculated the N fixation rates using the fodder crop production and fixation dataset of Einarsson et al.²¹ available at country-level during 1961–2019. The N fixation for the fodder crops in Einarsson et al.²¹ was estimated assuming a linear relationship between the crop yield and N fixed by the crops following the approach of Lassaletta et al.⁴⁷. We first derived the N fixation rates ($BNF_{fodder-Rate}$ ($kg \ kg^{-1}$ of crop product)) for the given fodder crop by dividing the values for the N fixation of fodder crops ($BNF_{fodder}$ ($kg \ yr^{-1}$)) with their respective production ($Pro_{fodder}$ ($kg \ yr^{-1}$)) as given in Eq. (67). The derived N fixation rates were then multiplied with the fodder crop production during 1850–2019 to get the N fixation by fodder crops ($BNF_{fodder}$ ($kg \ yr^{-1}$)) as expressed in Eq. (68).

$$BNF_{fodder-Rate}(c) = \frac{BNF_{fodder}(c)}{Pro_{fodder}(c)} \quad (67)$$

$$BNF_{fodder}(i, c, y_{1850-2019}) = Pro_{fodder}(i, c, y_{1850-2019}) \times BNF_{fodder-Rate}(c) \quad (68)$$

N fixation rates for pastures was taken as 5 kg of N $ha^{-1}$ as suggested by Bouwman et al.⁸. Moreover, we applied a N fixation rate of 1.77 kg $ha^{-1} \ yr^{-1}$ for boreal forest (for the forests in Norway, Sweden, Finland, Lithuania, Latvia, Estonia and Russia) and of 16 kg $ha^{-1} \ yr^{-1}$ for forest in temperate zone (for the forests in all other countries) as given by Cleveland et al.⁵⁴. For semi-natural vegetation, we applied the same N fixation rate given for the natural grassland following Cleveland et al.⁵⁴ (2.7 kg of N $ha^{-1} \ yr^{-1}$). N fixation rate for different crops and land types used in this study are listed in Table 2.

**N outputs.** The following section provides step wise reconstruction of the N removal from the croplands, pastures and forest. Furthermore, we provide details on the approach used for estimating gridded crop production which is utilized to derive the N removal from harvested crops.

**N removal from cropland.** N removal from cropland ($Rem_{cr}$ ($kg \ yr^{-1}$)) was derived as a sum of N removal from all crops calculated by multiplying the crop production ($Pro_{crops}$ ($t \ yr^{-1}$)) with crop specific N content ($N_{content}$ ($c$) ($kg \ t^{-1}$)); as given in Eq. (69). The coefficients applied to estimate N removal from harvested crops are given in the Table 2.

$$Rem_{cr}(i, y_{1850-2019}) = \sum_{i=1}^{n_c} (Pro_{crops}(i, c, y_{1850-2019}) \times N_{content}(c)) \quad (69)$$

Below, we further explained the steps to derive the crop specific annual production during the period 1850–2019.

**Crop production.** We obtained country-level datasets crop production from FAOSTAT for the period 1961–2019⁵⁶. FAOSTAT provides production dataset for 173 crops and excludes dataset for fodder crops. For fodder crops, we derived the dataset from Einarsson et al.²¹ that provides estimates for six fodder crops during 1961–2019 at country-level for European countries. For the missing countries and time period in dataset from Einarsson et al.²¹, we adapted the same methodology as explained in section “Crop specific harvested area-Fodder crops”. Prior to 1961, we reconstructed the country-level dataset from Bayliss-Smith and Wannali
(1984)⁵⁷ as cited in OWD⁵⁸ which provides wheat yield for the years 1850, 1909, 1934 and 1961 across a few European countries. We generated the annual time series of the wheat yield from OWD (denoted as Yield<sub>wheat,OWD</sub> (kg ha<sup>-1</sup> yr<sup>-1</sup>)) for the period 1850–1960 using linear interpolation. The time series of the wheat yield of the countries not included in the OWD dataset was computed as the mean of the wheat yield across the European countries for which data were available. After getting annual wheat yield from 1850 to 1960, we derived wheat production (Pro<sub>wheat</sub> (kg yr<sup>-1</sup>)) which is the product of the wheat yield (Yield<sub>wheat,OWD</sub> (kg ha<sup>-1</sup> yr<sup>-1</sup>)) and wheat harvested area (C<sub>A<sub>crops</sub></sub> (ha), calculated from Eq. (40).

For the reconstruction of crop production for all other crops (other than wheat) during the time period 1850–1960, we utilized the temporal variability of wheat production. To do so, we multiplied the country-level ratio between the wheat production from OWD (Pro<sub>wheat,OWD</sub>) for time period 1850–1960 and the wheat production from FAOSTAT (Pro<sub>wheat,FAO</sub>) in the year 1961 by the crop production (Pro) of all other crops derived from FAOSTAT (for non-fodder crops) and Einarsson et al.²¹ (for fodder crops) in 1961 as given in Eq. (70) (all units are in kg yr<sup>-1</sup>).

$$Pro_{\text{crops}}(u, c, y_{1850-1960}) = \frac{Pro_{\text{wheat,OWD}}(u, c, y_{1850-1960})}{Pro_{\text{wheat,FAO}}(u, c, y_{1961})} \times Pro(u, c, y_{1961}) \quad (70)$$

We downscaled the country-level crop productions (Pro<sub>crops</sub>) using the gridded database of crop production from Monfreda et al.³⁵ (Pro<sub>crops,Monfreda</sub>), as given in Eq. (71) (all units are in kg yr<sup>-1</sup>). Monfreda et al.³⁵ provides crop production for 175 different crops, including both fodder and non-fodder crops, globally at 5 arcmin spatial resolution for the year around 2000. As depicted in Eq. (71) below, this normalisation approach maintained the spatial heterogeneity and consistency in crop production estimates while taking into their respective annual trajectories at a country-level.

$$Pro_{\text{crops}}(i, c, y_{1850-2019}) = \frac{Pro_{\text{crops,Monfreda}}(i, c, y_{2000})}{\sum_{i=1}^{n_u} Pro_{\text{crops,Monfreda}}(i, c, y_{2000})} \times Pro_{\text{crops}}(u, c, y_{1850-2019}) \quad (71)$$

Here, u refers to a given country and n_u the total number of grid cells within a country.

**N removal from pasture.** We estimated N removal from pastures (Rem<sub>past</sub> (kg ha<sup>-1</sup> yr<sup>-1</sup>)) as harvested and grazed grass following a method proposed in previous studies⁸,⁵⁹. The method is chiefly based on a nitrogen use efficiency (NUE) in which N removal from pastures is derived by multiplying the N removal coefficient C<sub>Rem<sub>past</sub></sub> (-) with the N inputs to pastures (Inp<sub>past</sub> (kg ha<sup>-1</sup> yr<sup>-1</sup>)) adjusted for the losses of N after manure application to the soil (N<sub>losses</sub> (%)) as given in Eq. (72):

$$Rem_{\text{past}}(i, y_{1850-2019}) = C_{Rem_{\text{past}}} \times (Inp_{\text{past}}(i, y_{1850-2019}) - (N_{\text{losses}} \times Man_{\text{past}}(i, y_{1850-2019}))) \quad (72)$$

The N removal coefficient and N losses for the animal manure are varying among studies⁸,⁵⁹, and thus are highly uncertain. To account for these uncertainties, we considered two scenarios. In the first scenario, we assumed a value of 0.6 for C<sub>Rem<sub>past</sub></sub> based on Bouwman et al.⁸ and a value of 0.2 for N<sub>losses</sub> based on FAOSTAT. In the second scenario, we assumed a C<sub>Rem<sub>past</sub></sub> value of 0.4 and 0.5 for the countries located in Eastern and Western Europe, respectively, based on NUE values reported by Kaltenegger et al.⁶⁰ and 0 for N<sub>losses</sub>.

**N removal from forest.** N removal from forest (Rem<sub>for</sub> (kg ha<sup>-1</sup> yr<sup>-1</sup>)) was estimated by adapting the method from Chang et al.⁶¹, which assumes that N atmospheric deposition has a N fertilization effect and leads to an increase in woody biomass. To account for this aspect, we multiplied the atmospheric N deposition (DEP<sub>for</sub> (kg ha<sup>-1</sup> yr<sup>-1</sup>)) by the rate of forest N removal (C<sub>Rem<sub>for</sub></sub> (-)) as given in Eq. (73).

$$Rem_{\text{For}}(i, y_{1850-2019}) = C_{Rem_{\text{For}}} \times DEP_{\text{For}}(i, y_{1850-2019}) \quad (73)$$

We selected a value of 0.02 for C<sub>Rem<sub>For</sub></sub> based on the values suggested by Chang et al.⁶¹.

**N surplus across Europe (1850–2019).** We constructed 16 N surplus estimates at the spatial resolution of 5 arcmin (1/12°) across Europe, while accounting for the uncertainties in major components of N surplus. Specifically, we obtained two datasets for fertilizer, four datasets for animal manure and two datasets for N removal from pastures. Here, we present an average value of the 16 N surplus estimates during the time period of 1850–2019 in Figs. 2 and 3 and we show the uncertainties in Fig. 4.

The spatio-temporal variability in our total N surplus at gridded level is presented in Fig. 2. Here, we present the snapshot for the year 1900, 1930, 1960, 1990, 2000 and 2015 for the illustration. We observed a clear latitudinal pattern with lower values of the N surplus in northern Europe, higher values at mid-latitudes in particular in Western and Central Europe, to moderate value in southern Europe. Furthermore, N surplus showed large temporal changes during the period of 1850–2019 in most places, apart from northern Europe where its value remained relatively low and stable. This reflects the importance of having a long-term dataset of N surplus that allows quantification of these temporal developments. Up to 1900, N surplus in the majority of the grid cells did
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Snapshots of N surplus (kg ha⁻¹ of grid area yr⁻¹) across Europe. The figure shows the annual spatial variation in N surplus given as the mean of our 16 N surplus estimates for the selected years.</figcaption></figure>

not exceed 20 kg ha⁻¹ yr⁻¹. N surplus increased by a factor of 2–3 in around 70% of the grid cells during the time period of 1930–1990 (from around 20 kg ha⁻¹ yr⁻¹ in 1930 to values in the range 40–60 kg ha⁻¹ yr⁻¹ in 1990). In most of the grid cells of the industrialized European countries, N surplus reached its higher value around 1990 and declined in the subsequent time period. This decline was particularly apparent in countries such as Denmark, Germany, Latvia, Lithuanian, Estonia, Belarus, Ukraine etc. Conversely, Spain continued to show an upward trend in the majority of the grid cells during 1960–2000 with a gradual downwards trend afterwards. The gridded N surplus allows to examine the N surplus at sub-country levels, such as the levels of the European “Nomenclature of Territorial units for statistics” (NUTS) classification as depicted in Fig. 3. This figure shows the value of having information at sub-country level, since country-level (NUTS 0) averaging of N surplus can conceal large disparities between the territories. For instance, in 2015, the N surplus at NUTS 0 (country) level for France and the United Kingdom took moderate values (29 and 38 kg ha⁻¹ yr⁻¹, respectively). However, the NUTS 1 and NUTS 2 (sub-country) level information show that the N surplus was high in part of these countries, with values reaching 76 kg ha⁻¹ yr⁻¹ in the Brittany region in north-west France and 69 kg ha⁻¹ yr⁻¹ in Wales in the UK. Our gridded dataset also allows with the flexibility to assess the temporal trajectories of the N surplus for any river basin that may span different countries, as shown in Fig. 3 (right most panel). This is critical because river basins are relevant spatial scale for water quality modelling and land and water management¹²,⁶². Between 1930 and 1960, N surplus in the majority of the river basins showed a constant increasing pattern, a sharp jump in 1990 and declined afterwards. For instance, in the basin of the Danube river, that flows through much of the Central and Southeastern European countries, N surplus showed around a two-fold growth during 1930–1960 (14–26 kg ha⁻¹ yr⁻¹), further increased by a factor of 1.9 between 1960–1990 (49 kg ha⁻¹ yr⁻¹ in 1990), and reduced by around 39% during 1990–2015 (30 kg ha⁻¹ yr⁻¹ in 2015).

Figure 4 shows that the uncertainty range (between the maximum and minimum values of our 16 N surplus estimates) of agricultural and total N surplus estimates varies over time and across regions. Here for the illustration, we take N-surplus estimates for the scale of EU-28, Germany and Danube river basin. In general, for the agricultural N surplus, we observed that, during 1850–1910, the width of the uncertainty intervals (grey ribbons in Fig. 4a–c) has a similar magnitude to the mean estimate (red lines) for EU-28 and the Danube river and it is even more than twice the mean estimate for some years for Germany. Between 1910–2019, the width of the uncertainty range relative to the mean estimates decreased by a factor of around 7 for both EU-28 and the Danube river basin, i.e. it is equal to around 10–20% of the mean in 2019. For Germany, it took higher values than EU-28 and the Danube river basin for the recent years (around 40% of the mean estimate in 2019). In terms of absolute differences between the maximum and minimum estimates, we observed that the width of the uncertainty interval (grey ribbons) gradually increased during 1850–1950 (3–9 kg ha⁻¹ yr⁻¹) for EU-28, Germany and the Danube river basin (Fig. 4a–c). Between 1950–1990, the values for the Danube river basin kept on increasing (up to almost 15 kg ha⁻¹ yr⁻¹ in the 1980s), EU-28 showed a more stable pattern (10–13 kg ha⁻¹ yr⁻¹). After 1990, the uncertainty bound (ribbons) exhibited first a gradual decrease in thickness for EU-28 and the Danube river basin (down to around 5 kg ha⁻¹ yr⁻¹ in the 2000s) and a slight increase was seen in the recent years for EU-28. Comparatively, N surplus for Germany has a much larger uncertainty interval during 1950–2019 that keeps widening in time to reach a width of around 28 kg ha⁻¹ yr⁻¹ after 2015. We noticed a similar pattern in the uncertainty range over time for the total N surplus for EU28, Germany and the Danube river basin (Fig. 4d–f).
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Total N surplus (kg ha⁻¹ of physical area yr⁻¹) at multiple spatial levels for four years (1930, 1960, 1990, 2015). N surplus is given as the mean of our 16 N surplus estimates. NUTS: Nomenclature of Territorial units for statistics.</figcaption></figure>

These analyses highlight the value of our ensemble N surplus dataset, emphasising the uncertainty in N surplus estimates.

## Data Records

The dataset comprises 16 N surplus estimates in the Network Common Data Form (NetCDF) file format. Each NetCDF file contains an annual variable representing N surplus with the unit kg ha⁻¹ of grid area yr⁻¹ at 5 arcmin (1/12°) resolution for the period 1850–2019. We also provide N surplus estimates at aggregated spatial levels (i.e. NUTS 1, NUTS 2 and selected river basins covering in and around the areas of Europe) in comma-separated value (.csv) formats. The datasets are archived at the Zenodo repository⁶³. The description of the files provided in data repository is given below:

* 16 files termed as "N sur total kg ha grid 1850 2019 method xx v1.0" for gridded N surplus data in NetCDF format. Here, xx presents dataset from method 1 to method 16 and each of them contains 170 years (1850–2019) of gridded N surplus.

* 3 files termed as "N sur total kg ha NUTS xx" for N surplus at spatial aggregated levels in .csv format. Here, xx presents NUTS 1, 2 and 3 and each of them contains mean and standard deviation of 16 N surplus estimates.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Agricultural N surplus (kg ha⁻¹ of agricultural area yr⁻¹) and total N surplus (kg ha⁻¹ of physical area yr⁻¹) for EU-28, Germany and the Danube river basin (5 years moving average during 1850–2019). The grey color ribbon in each panel shows the ranges (minimum and maximum values) of the 16 N surplus estimates reconstructed in this study, whereas the average value is presented by a red line. (a–c) Agricultural N surplus for EU-28, Germany and Danube river, (d–f) Total N surplus for EU-28, Germany and Danube river.</figcaption></figure>

<table><caption>**Table 3.** Existing datasets used for technical evaluation and comparison in this study.</caption><thead><tr><td>Datasets</td><td>Budget type</td><td>Land type</td><td>Time period</td><td>Spatial information</td><td>References</td></tr></thead><tbody><tr><td>Eurostat</td><td>N land budget</td><td>Agriculture</td><td>1990–2019 (kg ha<sup>-1</sup> yr<sup>-1</sup> of agriculture area), 2010–2014 (tonne)</td><td>Country (Europe)</td><td>19</td></tr><tr><td>Zhang et al.<sup>5</sup></td><td>N soil surface budget</td><td>Cropland</td><td>1961–2015</td><td>Country (Global)</td><td>5</td></tr><tr><td>Hutchings et al.<sup>22</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>1990–2010</td><td>Country (DK)</td><td>22</td></tr><tr><td>Defra<sup>23</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>1990, 1995 and 2000 to 2020</td><td>Country (UK)</td><td>23</td></tr><tr><td>Poisvert et al.<sup>24</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>1940–2010</td><td>Country (FR)</td><td>24</td></tr><tr><td>Eriksson et al.<sup>25</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>1960–2000</td><td>Country (PO)</td><td>25</td></tr><tr><td>Häussermann et al.<sup>26</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>1995–2017</td><td>Country (DE)</td><td>26</td></tr><tr><td>Behrendt et al.<sup>64</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>1950–1998</td><td>NUTS 1 (DE)</td><td>64</td></tr><tr><td>Leip et al.<sup>20</sup></td><td>N soil surface budget</td><td>Agriculture</td><td>2001–2003</td><td>Country (Europe)</td><td>20</td></tr></tbody></table>

* 1 file termed as “N_sur_total_kg_ha_river_basin” for N surplus at selected river basin levels in .csv format. This file contains mean and standard deviation of 16 N surplus estimates.

Additionally, readme files are provided for the NUTS and river basins ID's.

## Technical Evaluation

To check the spatio-temporal plausibility and consistency of our N surplus dataset, we compared our 16 N surplus estimates with available datasets (see Table 3), due to a lack of direct observations of N surplus. First, we used the datasets of the statistical office of the European Union (Eurostat)¹⁹ and Zhang et al.⁵, available at country-level across Europe for the period 1990–2019 and 1961–2015, respectively. The dataset of Eurostat¹⁹ is a reference dataset for Europe that consistently provides N land budgets for agricultural areas across 37 European countries. Zhang et al.⁵ provides different components of the N surplus for croplands by synthesising information of 13 different N datasets; and as such it has been proposed as a benchmark estimate⁵. To check for further consistency, we compared our N-budget with country specific estimates for Denmark²² (1990–2010), the United Kingdom²³ (1990, 1995 and 2000–2018), France²⁴ (1940–2010), Poland²⁵ (1960–2000) and Germany²⁶,⁶⁴, (1950–2017), while results for 25 European countries were additionally compared to estimates given by Leip et al.²⁰ for the period of 2001–2003. For most of these studies, we could only extract information at country-level. An exception is the German dataset for which we could derive NUTS 1 (state) level data for the period 1950–2015 combining the two datasets of Behrendt et al.⁶⁴ and Häussermann et al.²⁶ (for methodological details we refer to Ebeling et al.⁶⁵).
At continental-level, the uncertainty intervals of our N surplus (assessed as one standard deviation around the mean of the 16 estimates) captured the values derived from Eurostat¹⁹ and Zhang et al.⁵. Specifically, our average cropland N surplus for 20 countries over the period 1961–2015 is equal to 7.86 ± 1.18 Tg of N, whereas Zhang et al.⁵ provides a value of 9.2 Tg of N. Our average agricultural N budget for EU-28 over the period 2010–2014 is equal to 8.9 ± 0.75 Tg of N, while the value calculated from Eurostat is 8.88 Tg of N. Moreover, with respect to the Leip et al.²⁰ study which reported an average N surplus of 59 kg ha⁻¹ yr⁻¹ for the period 2001–2003 for EU-27, our dataset show a reasonably good match with N surplus values of around 60 ± 4 kg ha⁻¹ yr⁻¹ for the same time period. With respect to country specific study, Häussermann et al.²⁶ reported an average value of 77 kg ha⁻¹ yr⁻¹ for N soil surface budget for Germany during 2015–2017, which is also included in our uncertainty bound of 70 ± 9 kg ha⁻¹ yr⁻¹. N surplus estimates for both France (47 ± 4 kg ha⁻¹ yr⁻¹, average over 1940–2010) and Poland (55 ± 4 kg ha⁻¹ yr⁻¹ an average of 1960–2000) are slightly higher in our study compared to the values found in previous studies that are equal to 37 kg ha⁻¹ yr⁻¹²⁴ and 41 kg ha⁻¹ yr⁻¹²⁵, respectively.

To determine the differences in the average N surplus from this study with respect to Eurostat data (respectively Zhang et al.⁵), we calculated the country-wise relative difference *Diff* (%), as given in Eq. 74:

$$Diff = \left( \frac{Reference_{Nsur} - Study_{Nsur}}{Reference_{Nsur}} \right) \times 100 \qquad (74)$$

where *Reference*<sub>Nsur</sub> denotes the average N budget given by Eurostat during the time period 1990–2019 (in kg ha⁻¹ of agricultural area yr⁻¹) (respectively Zhang et al.⁵ between 1961–2015 in tonne) and *Study*<sub>Nsur</sub> denotes one of our 16 N surplus estimates averaged over the relevant time period. We present here the mean and standard deviation (SD) for the 16 value of *Diff* calculated for each of our 16 N surplus datasets. We found that, in about 70% of the countries, the mean *Diff* between our estimates and Eurostat was contained in the range ± 30% (see Fig. 5c). Large values of the mean *Diff* (exceeding ± 50%) were observed for some countries, namely Romania (RO; *Diff* = −111), Ireland (IE; *Diff* = −63) and Greece (EL; *Diff* = 50). We note that these differences in N surplus between Eurostat and our estimates could occur because Eurostat unlike our estimate accounts for volatility losses during manure management and storage. With respect to Zhang et al.⁵, 85% of the countries show a mean *Diff* in the range ± 30% (see Fig. 6c), and all countries have a mean *Diff* within the range ± 50%. In particular, higher values (in the range 32–50%) are found for the United Kingdom (UK), France (FR) and Austria (AT). Importantly, when considering the entire interval of *Diff* (across our 16 estimates) beyond the mean value only, we found a better correspondence between Eurostat¹⁹ and Zhang et al.⁵ datasets and ours. In fact, we observed that the interval contains values in the range ± 30% for 80% of the countries with respect to Eurostat¹⁹ data and for 95% of countries with respect to Zhang et al.⁵ data. Nonetheless, in general, the bounds of our 16 values of *Diff* are reasonably narrow for both Eurostat¹⁹ and Zhang et al.⁵. In particular, we observed that the SD of *Diff* is lower than or equal to around 11% for all countries (except for Romania) with respect to Eurostat¹⁹ and for 70% of the countries with respect to Zhang et al.⁵ data.

To assess the consistency of the temporal dynamics between the N budget from this study and available datasets, we computed a correlation coefficient (*r*). Here, we present the mean and standard deviation of 16 *r* values derived for each of our 16 N surplus datasets. Overall, we found a strong, positive correlation between estimates of this study and Eurostat¹⁹ during 1990–2015 (*r* = 0.90 ± 0.03; see Fig. 5b). In three-quarters of the countries, we observed high mean *r* values in the range 0.6–0.9, and none of the countries have a mean *r* value below 0.5 (Fig. 5a). The dynamics of our cropland N budgets also shows positive and strong correlation with Zhang et al.⁵ during 1961–2015 (*r* = 0.95 ± 0.02); see Fig. 6b). In around three-quarters of the countries, there is a very good consistency between our estimates and Zhang et al.⁵, with mean *r* values in the range 0.8–0.99 (Fig. 6a). In the Netherlands (NL), Finland (FI) and the United Kingdom (UK), we noticed lower mean values of *r* (in the range 0.3–0.6). However, the SD of the *r* values across the 16 N surplus estimates is large for these three countries. This results in values of the upper bound of *r* (one SD above the mean value) equal to 0.75, 0.56 and 0.77 for NL, FI and UK, respectively, which is much higher than the mean *r* value. This means, that the temporal dynamics of some of our N surplus estimates are reasonably consistent with Zhang et al.⁵ data for these three countries. Overall for most countries, we found that the spread of our 16 *r* values is reasonably narrow for both Eurostat¹⁹ and Zhang et al.⁵. This reflects that the temporal dynamics across our 16 N surplus estimates are similar. Specifically, the SD of *r* is lower than 0.10 for 80% of the countries with respect to Eurostat¹⁹ and for 85% of the countries with respect to Zhang et al.⁵.

Similarly, we observed that our N budget is strongly and positively correlated with estimates for EU countries (*r* = 0.94 ± 0.01) during 2001–2003 compared to Leip et al.²⁰ estimates (see Fig. 7a), Denmark (DK) (*r* = 0.98 ± 0.01) between 1990–2010 (see Fig. 7b), and United Kingdom (UK) (*r* = 0.96 ± 0.01) for the time period 1990–2009 (Fig. 7c) and Germany (DE) (*r* = 0.96 ± 0.01) between 1950–2015 (Fig. 7d). Furthermore, for Germany we compared our estimates at a sub-country (NUTS 1) level based on the harmonized dataset of Behrendt et al.⁶⁴ and Häussermann et al.²⁶ for the years 1950–2015. We examined the results aggregated over the states of the former West and East Germany and found *r* values of 0.93 ± 0.04 and 0.88 ± 0.04, respectively, showing a strong and positive association between the sub-country scale German dataset and ours (see Fig. 7e,f). We also observe for some years much lower values of the N surplus for East Germany in the dataset²⁶,⁶⁴ compared to ours. The two datasets show different temporal dynamics around the year of the German reunification (1990), when agricultural activities were disrupted.

Overall, we found that our 16 estimates of N surplus at continental and country-level captured reasonably well other existing estimates provided by previous studies in terms of their spatio-temporal pattern and relative differences. However, we found some differences in N surplus estimates, which emphasizes the need for a more
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Country level comparison between N surplus for agriculture soils estimated by Eurostat¹⁹ and this study for the period 1990–2019. The circles in panels (b) and (c) denote the average of the 16 N surplus estimates reconstructed in this study, whereas the bars show the standard deviation. (a) Pearson correlations coefficients (r) values for every country. Countries with white color are excluded from the comparison because they are not part of the Eurostat dataset. (b) Linear fit between the N surplus values for all countries and all years in the two studies: x-axis shows the N surplus calculated in this study and y-axis presents the N surplus given by Eurostat, (c) Relative difference (defined in Eq. 74) in N surplus in this study with respect to Eurostat against correlation coefficient for each country. Romania (RO) is not presented in this plot as it is an outlier with a value of −111% for mean relative difference and a SD of 35% around the mean. The correlation for Romania is however high (0.9 ± 0.03).</figcaption></figure>

comprehensive characterization of the uncertainties in the N surplus in the future, to further improve the N surplus estimation. In fact, N surplus datasets are the result of a modeling process involving the use of uncertain underlying data and numerous uncertain methodological choices. Only part of these uncertainties were considered in our study, while we discuss in the following other sources of uncertainty that could explain the differences between our estimates and previous studies. First, the input data sources of different components of the N surplus such as fertilizer, animal manure, atmospheric deposition differ among studies, which could lead to differences in estimated amounts among studies. For example, regarding fertilizers, the database of IFA⁶⁶ relies on sales/trades data that mainly come from fertilizer companies of their member countries, FAOSTAT⁴¹ collects information officially from its member countries via questionnaire and Eurostat¹⁹ obtained consumption amount from country specific statistics (i.e. data based on trade, production of fertilizer and surveys) and from UNFCCC (United Nations Framework Convention on Climate Change) for the missing countries⁶⁷. Second, studies made various assumptions regarding the amount of fertilizer and animal manure applied to cropland and pastures. For example, many studies documented in Zhang et al.⁵ did not exclude the fertilizer used on pastures to assess cropland application. Thus, there is a large uncertainty in the spatial distribution of applied N. Moreover, differences in N budget can largely be attributed to different coefficients used by different studies for estimating for instance N content in harvested crops, N fixation and volatilization losses in animal manure²⁸. Another plausible cause of uncertainty in N budget is the different crop types considered in the studies, especially fodder crops²⁷.

Despite the above-mentioned challenges, our dataset is a step forward since it provides 16 uncertainty estimates in major components of the N surplus, while previous datasets⁵,¹⁹,²¹⁻²³,²⁵,²⁶ typically do not report uncertainties in N surplus. The findings of the comparison of our N surplus estimates with previous datasets highlight the need to further characterize the uncertainties in future versions of the N surplus dataset.
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Country level comparison between N surplus for the cropland estimated by Zhang *et al.*⁵ and this study for the period 1961–2015. The circles in panel (b) and (c) denote the average of the 16 N surplus estimates reconstructed in this study, whereas the bars show the standard deviation. (a) Pearson correlations coefficients (r) values for every country. Countries with white color are excluded from the comparison because they are not part of the Zhang *et al.*⁵ dataset, (b) Linear fit between the N surplus values for all countries and all years in the two studies: x-axis shows the N surplus calculated in this study and y-axis presents the N surplus given by Zhang *et al.*⁵(c) Relative difference (defined in Eq. 74) in N surplus in this study with respect to Zhang *et al.*⁵ against correlation coefficient for each country.</figcaption></figure>

## Usage Notes

This study makes available annual estimates of N surplus from 1850 to 2019 at gridded resolution across Europe. We accounted for uncertainties due to different underlying data and methodological choices in components of N surplus that are reported to have considerable uncertainties³, namely N inputs from fertilizer and manure, and N removal from pastures. Our uncertainty ensemble was able to generally capture well country-level N surplus estimates available from previous studies/reports/databases (see Section “Technical evaluation”). We do acknowledge that our dataset also does not include all sources of uncertainties, and we discussed further uncertain data and choices that are worth investigating in future versions of N surplus estimates (Section “Technical evaluation”). In particular, further considerations can be given to better inform about disaggregation schemes and their associated uncertainties. In this study, we downscaled the country level estimates (e.g. the fertilizer application amounts for the different crop categories in 2014/2015 from IFA⁴⁰) using the spatial variability in land use areas (Ramankutty *et al.*³⁴ dataset), and crop specific harvested areas and crop production from the Monfreda *et al.*³⁵ dataset. Moreover, the construction of country-level N surplus components before 1961 is based on either global temporal dynamics (e.g. for fertilizers) or proxy information (e.g. for crop yields, we assumed the same dynamic of yield as of wheat). Hence, our database is likely to have less reliability during the time period 1850–1960. Comparatively, we can consider that our estimates have a lower uncertainty and a moderate reliability after 1961, since there is better availability of N surplus components at country level (e.g. estimates of fertilizer and manure from FAOSTAT). Furthermore, we expect that our N surplus data are less uncertain for the recent time period (from the mid-1990s). This is mainly due to the fact that we utilized detailed spatial dataset for the reconstruction of N surplus estimates during this time period (e.g. land use and crop production datasets that are available around 2000).

Therefore, users should be aware that our N surplus data are more reliable for the recent time frame, while past estimates are valuable for informing general trends in the N surplus values. Regarding the spatial scale, we recommend not to use our dataset at gridded level, because of the uncertainties in the spatial disaggregation schemes, but rather at higher aggregation levels, such as country level, the different European socio-economic regions/states/provinces (e.g. NUTS 0, NUTS 1, NUTS 2 levels), and river basin scale (see Section “N budgets across Europe” and Fig. 3) to support different land and water management activities. Importantly, we provide
<figure><img src="image_9.png" /><figcaption>Fig. 7 Scatter plots between the N surplus for agricultural soil from various previous studies and our study for different countries and corresponding linear fits. The circles in each panel denote the average of 16 N surplus reconstructed in this study, whereas the bars show the standard deviation. (a) European countries during 2001–2003 (each point in this plot represent one country), (b) Denmark for each year in 1990–2010, (c) United Kingdom (UK) for 1990, 1995 and each year in 2000–2018, (d) Germany (DE) for each year in 1950–2015, (e) Aggregated values over the states of the former West Germany (DE-West) for each year in 1950–2015, (f) Aggregated values over the states of the former East Germany (DE-East) for each year in 1950–2015.</figcaption></figure>

here a consistent methodology using publicly available information for the reconstruction of N surplus over a long time period. Our methodology thus, can be reused to incorporate new underlying datasets and further updates of datasets that we already used, such as a new version of the FAOSTAT database, as they become available.

## Code availability

We used the RStudio version (1.2.5019) for data processing. The codes for the study are available at GitLab https://git.ufz.de/batool/n_s surplus.

Received: 28 May 2022; Accepted: 4 September 2022;

Published online: 10 October 2022

## References

1. Lu, C. & Tian, H. Global nitrogen and phosphorus fertilizer use for agriculture production in the past half century: Shifted hot spots and nutrient imbalance. *Earth Syst. Sci. Data* **9**, 181–192, https://doi.org/10.5194/essd-9-181-2017 (2017).

2. Bodirsky, B. L. et al. Reactive nitrogen requirements to feed the world in 2050 and potential to mitigate nitrogen pollution. *Nat. Commun.* **5**, https://doi.org/10.1038/ncomms4858 (2014).

3. Dentener, F. et al. Nitrogen and sulfur deposition on regional and global scales: A multimodel evaluation. *Global Biogeochem. Cycles* **20**, https://doi.org/10.1029/2005GB002672 (2006).

4. Galloway, J. N. et al. Galloway-et-al-2004-BioGeoCh, vol. 70 (2004).

5. Zhang, X. et al. Quantification of global and national nitrogen budgets for crop production. *Nat. Food* **2**, 529–540, https://doi.org/10.1038/s43016-021-00318-5 (2021).

6. Gobler, C. J. Climate Change and Harmful Algal Blooms: Insights and perspective. *Harmful Algae* **91**, 101731, https://doi.org/10.1016/j.hal.2019.101731 (2020).

7. Erisman, J. W. et al. Consequences of human modification of the global nitrogen cycle. *Philos. Trans. R. Soc. B Biol. Sci.* **368**, https://doi.org/10.1098/rstb.2013.0116 (2013).

8. Bouwman, A. F., Van Brecht, G. & Van Der Hoek, K. W. Global and regional surface nitrogen balances in intensive agricultural production systems for the period 1970–2030. *Pedosphere* **15**, 137–155 (2005).

9. WHO. Nitrate and nitrite in drinking water: background document for development of WHO guidelines for drinking water quality. *Geneva World Heal. Organ.* **31** (2016).
10. Wang, L., Butcher, A. S., Stuart, M. E., Gooddy, D. C. & Bloomfield, J. P. The nitrate time bomb: A numerical way to investigate nitrate storage and lag time in the unsaturated zone. *Environ. Geochem. Health* **35**, 667–681, https://doi.org/10.1007/s10653-013-9550-y (2013).

11. Eurostat. *Agriculture, forestry and fishery statistics: 2020 edition* (2020).

12. Basu, N. B. et al. Managing nitrogen legacies to accelerate water quality improvement. *Nat. Geosci.* **15**, 97–105, https://doi.org/10.1038/s41561-021-00889-9 (2022).

13. Ascott, M. J. et al. The need to integrate legacy nitrogen storage dynamics and time lags into policy and practice. *Sci. Total Environ.* **781**, 146698, https://doi.org/10.1016/j.scitotenv.2021.146698 (2021).

14. Vero, S. E. et al. Revue: L'état environnemental et les implications du décalage temporel du nitrate en Europe et Amérique du Nord. *Hydrogeol. J.* **26**, 7–22, https://doi.org/10.1007/s10040-017-1650-9 (2018).

15. McDowell, R. W., Simpson, Z. P., Ausseil, A. G., Etheridge, Z. & Law, R. The implications of lag times between nitrate leaching losses and riverine loads for water quality policy. *Sci. Rep.* **11**, 1–14, https://doi.org/10.1038/s41598-021-95302-1 (2021).

16. Byrnes, D. K., Van Meter, K. J. & Basu, N. B. Long-Term Shifts in U.S. Nitrogen Sources and Sinks Revealed by the New TREND-Nitrogen Data Set (1930–2017). *Global Biogeochem. Cycles* **34**, 1–16, https://doi.org/10.1029/2020GB006626 (2020).

17. Oenema, O., Kros, H. & De Vries, W. Approaches and uncertainties in nutrient budgets: Implications for nutrient management and environmental policies. *Eur. J. Agron.* **20**, 3–16, https://doi.org/10.1016/S1161-0301(03)00067-4 (2003).

18. FAOSTAT (Food and Agriculture Organization Corporate Statistical Database): Soil nutrient budget, available at: https://www.fao.org/faostat/en/#data/ESB [Accessed: 20-September-2021].

19. Eurostat: Gross nutrient balance on agricultural land, available at: https://ec.europa.eu/eurostat/datab browser/product/page/T2020_RN310 [Accessed: 10-September-2021].

20. Leip, A., Britz, W., Weiss, F. & De Vries, W. Farm, land, and soil nitrogen budgets for agriculture in Europe calculated with CAPRI. *Environ. Pollut.* **159**, 3243–3253, https://doi.org/10.1016/j.envpol.2011.01.040 (2011).

21. Einarsson, R. et al. Crop production and nitrogen use in European cropland and grassland 1961–2019. *Scientific Data* **1–29**, https://doi.org/10.1038/s41597-021-01061-z (2021).

22. Hutchings, N. J. et al. A nitrogen budget for Denmark; Developments between 1990 and 2010, and prospects for the future. *Environ. Res. Lett.* **9**, https://doi.org/10.1088/1748-9326/9/11/115012 (2014).

23. Defra: Soil Nutrient Balances, available at: https://www.gov.uk/government/statistics/uk-and-england-soil-nutrient-balances-2020 [Accessed: 15-February-2022].

24. Poisvert, C., Curie, F. & Moatar, F. Annual agricultural N surplus in France over a 70-year period. *Nutr. Cycl. Agroecosystems* **107**, 63–78, https://doi.org/10.1007/s10705-016-9814-x (2017).

25. Eriksson, H., Pastuszak, M., Löfgren, S., Mörth, C. M. & Humborg, C. Nitrogen budgets of the Polish agriculture 1960-2000: Implications for riverine nitrogen loads to the Baltic Sea from transitional countries. *Biogeochemistry* **85**, 153–168, https://doi.org/10.1007/s10533-007-9126-y (2007).

26. Häußermann, U. et al. Nitrogen soil surface budgets for districts in Germany 1995 to 2017. *Environ. Sci. Eur.* **32**, https://doi.org/10.1186/s12302-020-00382-x (2020).

27. Bouwman, L. et al. Exploring global changes in nitrogen and phosphorus cycles in agriculture induced by livestock production over the 1900-2050 period. *Proc. Natl. Acad. Sci. USA* **110**, 20882–20887, https://doi.org/10.1073/pnas.1012878108 (2013).

28. Bodirsky, B. L. et al. N2O emissions from the global agricultural nitrogen cycle-current state and future scenarios. *Biogeosciences* **9**, 4169–4197, https://doi.org/10.5194/bg-9-4169-2012 (2012).

29. Zhang, B. et al. Manure nitrogen production and application in cropland and rangeland during 1860&ndash;2014: A 5-minute gridded global data set for Earth system modeling. *Earth Syst. Sci. Data Discuss.* **1–35**, https://doi.org/10.5194/essd-2017-11 (2017).

30. West, P. C. et al. Leverage points for improving global food security and the environment. *Science* (80–). **345**, 325–328, https://doi.org/10.1126/science.1246067 (2014).

31. Potter, P., Ramankutty, N., Bennett, E. M. & Donner, S. D. Characterizing the spatial patterns of global fertilizer application and manure production. *Earth Interact.* **14**, https://doi.org/10.1175/2009EI288.1 (2010).

32. Xu, R. et al. Increased nitrogen enrichment and shifted patterns in the world's grassland: 1860-2016. *Earth Syst. Sci. Data* **11**, 175–187, https://doi.org/10.5194/essd-11-175-2019 (2019).

33. FAOSTAT: Food and Agriculture Organization Corporate Statistical Database, available at: https://www.fao.org/faostat/en/#home.

34. Ramankutty, N., Evan, A. T., Monfreda, C. & Foley, J. A. Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000. *Global Biogeochem. Cycles* **22**, 1–19, https://doi.org/10.1029/2007GB002952 (2008).

35. Monfreda, C., Ramankutty, N. & Foley, J. A. Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. *Global Biogeochem. Cycles* **22**, 1–19, https://doi.org/10.1029/2007GB002947 (2008).

36. Goldewijk, K. K., Beusen, A., Doelman, J. & Stehfest, E. Anthropogenic land use estimates for the Holocene - HYDE 3.2. *Earth Syst. Sci. Data* **9**, 927–953, https://doi.org/10.5194/essd-9-927-2017 (2017).

37. Goldewijk, K. K. *Estimating global land use change over the past 300 years: The HYDE database* https://doi.org/10.1029/1999GB001232 (2001).

38. FAOSTAT (Food and Agriculture Organization Corporate Statistical Database): Land Use domain, available at: https://www.fao.org/faostat/en/#data/RL/ [Accessed: 10-August-2021].

39. Bartholomé, E. & Belward, A. S. GLC2000: A new approach to global land cover mapping from earth observation data. *Int. J. Remote Sens.* **26**, 1959–1977, https://doi.org/10.1080/01431160412331291297 (2005).

40. Heffer, P. Assessment of Fertilizer Use by Crop at the Global Level. *Int. Fertil. Ind. Assoc.* **5**, 9 (2013).

41. FAOSTAT (Food and Agriculture Organization Corporate Statistical Database): Fertilizer by nutrients, available at: https://www.fao.org/faostat/en/#data/RFN [Accessed: 10-October-2021].

42. Holland, E. A., Lee-Taylor, J., Nevison, C. & Sulzman, J. Global N Cycle: Fluxes and N 2 O Mixing Ratios Originating from Human Activity. *Report* **4** (2005).

43. Turner, B. et al. *The Earth as transformed by human action: global change and regional changes in the biosphere over the past 300 years* (Cambridge University Press, with Clark University, 1990).

44. Nevison, C. D., Esser, G. & Holland, E. A global model of changing n2o emissions from natural and perturbed soils. *Climatic Change* **32**, 327–378 (1996).

45. FAOSTAT (Food and Agriculture Organization Corporate Statistical Database): Livestock Manure, available at: https://www.fao.org/faostat/en/#data/EMN [Accessed: 10-October-2021].

46. Dong, H. et al. 2006 IPCC Guidelines for National Greenhouse Gas inventories - Chapter 10: Emissions from livestock and manure management. *Report* (2020).

47. Lassaletta, L., Billen, G., Grizzetti, B., Anglade, J. & Garnier, J. 50 year trends in nitrogen use efficiency of world cropping systems: The relationship between yield and nitrogen input to cropland. *Environ. Res. Lett.* **9**, https://doi.org/10.1088/1748-9326/9/10/105011 (2014).

48. FAO. *Nitrogen inputs to agricultural soils from livestock manure. New statistics*, vol. **24** (2018).

49. Robinson, T. P. et al. Mapping the global distribution of livestock. *PLoS One* **9**, https://doi.org/10.1371/journal.pone.0096084 (2014).

50. Bouwman, A. F. et al. Lessons from temporal and spatial patterns in global use of N and P fertilizer on cropland. *Sci. Rep.* **7**, 1–11, https://doi.org/10.1038/srep40366 (2017).
51. Tian, H. et al. The global N2O model intercomparison project. *Bull. Am. Meteorol. Soc.* **99**, 1231–1251, https://doi.org/10.1175/BAMS-D-17-0212.1 (2018).

52. Kremer, A. M. Methodology and Handbook Eurostat/OECD Nutrient Budgets, version 1.02. *Report* 1–112 (2013).

53. input4mips (input datasets for Model Intercomparison Projects): N deposition, available at: https://esgf-node.llnl.gov/search/input4mips/ [Accessed: 10-July-2021].

54. Cleveland, C. C. et al. Global patterns of terrestrial biological nitrogen (N2) fixation in natural ecosystems. *Global Biogeochem. Cycles* **13**, 623–645, https://doi.org/10.1002/(ISSN)1944-9224 (1999).

55. Van Meter, K. J., Basu, N. B. & Van Cappellen, P. Two centuries of nitrogen dynamics: Legacy sources and sinks in the Mississippi and Susquehanna River Basins. *Global Biogeochem. Cycles* **31**, 2–23, https://doi.org/10.1002/2016GB005498 (2017).

56. FAOSTAT (Food and Agriculture Organization Corporate Statistical Database): Crop production quantity, available at: https://www.fao.org/faostat/en/#data/QCL. [Accessed: 10-August-2021].

57. Bayliss-Smith & Wanmali. *Understanding Green Revolutions: Agrarian Change and Development Planning in South Asia* (Cambridge University Press, 1984).

58. Our World in Data: Crop Yields, available at: https://ourworldindata.org/crop-yields [Accessed: 10-April-2021].

59. Kaltenegger, K. & Winiwarter, W. Global Gridded Nitrogen Indicators: Influence of Crop Maps. *Global Biogeochem. Cycles* **34**, 1–17, https://doi.org/10.1029/2020GB006634 (2020).

60. Kaltenegger, K., Erb, K. H., Matej, S. & Winiwarter, W. Gridded soil surface nitrogen surplus on grazing and agricultural land: Impact of land use maps. *Environ. Res. Commun.* **3**, https://doi.org/10.1088/2515-7620/abedd8 (2021).

61. Chang, S. Y., Zhang, Q., Byrnes, D. K., Basu, N. B. & van Meter, K. J. Chesapeake legacies: The importance of legacy nitrogen to improving Chesapeake Bay water quality. *Environ. Res. Lett.* **16**, https://doi.org/10.1088/1748-9326/ac0d7b (2021).

62. Sarrazin, F. J. et al. Characterizing Catchment-Scale Nitrogen Legacies and Constraining their Uncertainties. *Water Resour. Res.* **1–32**, https://doi.org/10.1029/2021wr031587 (2022).

63. Batool, M., Sarrazin, F. J., Attinger, S. & Kumar, R. Long-term annual soil nitrogen surplus across europe (1850–2019). *Zenodo* https://doi.org/10.5281/zenodo.6581441 (2022).

64. Behrendt, H. et al. *Nutrient Emissions into River Basins of Germany on the Basis of a Harmonized Procedure (Research report 299 22 285).* Texte | 82/2003 Last access: 07 November 2021 (German Federal Environment Agency (Umweltbundesamt), Berlin, Germany, 2003).

65. Ebeling, P. et al. Water quality, discharge and catchment attributes for large-sample studies in Germany–quadica. *Earth System Science Data Discussions* **2022**, 1–41, https://doi.org/10.5194/essd-2022-6 (2022).

66. Mueller, N. D. et al. Closing yield gaps through nutrient and water management. *Nature* **490**, 254–257, https://doi.org/10.1038/nature11420 (2012).

67. Eurostat: Gross nutrient balance, available at: https://ec.europa.eu/eurostat/cache/metadata/en/aei_pr_gnb_esms.htm [Accessed: 10-October-2021].

68. Smil, V. Nitrogen in crop production: An account of global flows adds. *Global Biogeochem. Cycles* **13**, 647–662 (1999).

## Acknowledgements

The authors gratefully acknowledge the Global Water Quality Analysis and Service Platform (GlobeWQ) project for funding this work financed by the German Ministry for Education and Research (grant number 02WGR1527A). We are also thankful to UFZ for providing computing power and technical support to EVE super computing facility.

## Author contributions

R.K. and F.S. designed the study with inputs from N.B. and K.V.; and M.B. collected data and conducted the analysis with support from F.S. and R.K. M.B. drafted the initial manuscript together with F.S. F.S. and R.K. edited and revised the manuscript. All authors provided feedback and contributed to finalizing the manuscript.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to M.B. or R.K.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022