{'fraction_non_alphanumeric': 0.04357332500390442, 'fraction_numerical': 0.031183299495028374, 'mean_word_length': 4.586815317498788, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 27, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long-term annual soil nitrogen surplus across Europe (1850–2019)

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) 1,2 and deposition of atmospheric N originating from fossil fuel combustion [3][4][5] . 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 [6][7][8][9] . 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 10,11 . 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" 12,13 . These legacies are largely responsible for time lags between the implementation of water quality measures and the water quality response 14,15 .
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 16 . 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) 11 , or the soil surface budget which further excludes gaseous losses during manure management and storage 17 . 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 18 and for different underlying components of the agricultural N budget since 1961. More recently, Zhang et al. 5 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 19 ; Leip et al. 20 ) or to croplands only (Einarsson et al. 21 ). Moreover, efforts have been expanded to provide N budgets at sub-national levels for individual countries, such as, Denmark 22 , United Kingdom 23 , France 24 , Poland 25 and Germany 26 . Notably, uncertainty estimation is also currently not a common practice in N budget construction 17 , although uncertainties in existing N datasets are understood to be large 5 . 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 5,27,28 .
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 11,20,[22][23][24][25][26]29,30 . 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 1, 16,31,32 . 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. 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 11,20,26 . 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 −1 yr −1 ): soil agri NonAgri where i is grid cell; y is year; Surp soil is the total N surplus; Surp agri is the N surplus over agriculture area; and Surp NonAgri 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 cr ) and pastures (Surp past ). N surplus over croplands and pastures is estimated as the difference between the inputs to cropland and pasture (Inp cr and Inp past , respectively) through mineral fertilizers (FERT cr and FERT past , respectively), animal manure (MAN cr and MAN past , respectively), biological N fixation (BNF cr and BNF past , respectively) and atmospheric deposition (DEP cr and DEP past , respectively), and the N outputs (Out cr and Out past , respectively) through N removal via harvested crops (Rem cr ) and via animal grazing and cutting of grass (Rem past ), respectively. These can be summarized in Eqs. (2)(3)(4)(5)(6)(7)(8) (all variables are in kg ha −1 yr −1 ) as: past past past past past NatVeg NatVeg NatVeg NonVeg NonVeg Urban Urban = 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 11,34,38 . Pasture area in our estimates is the land under permanent meadow and pasture 38 . First, to obtain the spatial variability at gridded level, we used the dataset from Ramankutty et al. 34  www.nature.com/scientificdata www.nature.com/scientificdata/ each year during 1850-2019, by referencing it to the year 2000 (R HYDE-cr (-) and R HYDE-past (-), respectively), as given in Eqs. (16) and (17):  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. 34 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. 34 , while considering the annual temporal developments from HYDE. Finally, we derived agriculture area (A agri (ha)) by combining our reconstructed gridded cropland and pasture area (A cr (ha) and A past (ha), respectively) for the time period 1850-2019, as given in Eq. (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 cr (-) and R A past (-), respectively) between those given in FAOSTAT (A FAO-cr (ha) and A FAO-past (ha), respectively) and ones estimated in our study (A cr (ha) and A past (ha)), as summarized in Eqs. (21)(22):  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 ha ( ) , respectively) can be given in Eqs. (23)(24) as:  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):

= ×
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): www.nature.com/scientificdata www.nature.com/scientificdata/ ) 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. 34 ) 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 ha ( ) A cr ) and pasture areas (C ha ( ) A past ) 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 NonAgri (ha) that is left after allocation of cropland and pasture area, as given in Eq. (28):

= − +
Afterwards, we utilized the classification of land cover categories from global land cover (GLC) 39 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 (For GLC , NatVeg GLC , NonVeg GLC , Urban GLC , (-) respectively) from GLC land cover within our target grid cell of 5 arcmin, and their sum which corresponds to the total non-agricultural fraction (NonAgri GLC (-)). We then estimated the relative proportion of forest, semi-natural vegetation, non-vegetation and urban area fr For GLC , fr NatVeg GLC , fr NonVeg GLC , fr Urban GLC , (-) respectively) within the total non-agricultural fraction as given from Eqs. (29) to (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):     35 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) 35 . Among these 175 crops, we selected 17 major crops for which fertilizer application rates are provided 40 and that have either a large production across Europe or a N content greater than 1 kg of N tonne −1 of product. The selected crops are listed in Table 2.
www.nature.com/scientificdata www.nature.com/scientificdata/ To construct annual estimates of crop specific harvested area (A crops (ha)), we first applied temporal dynamics of cropland area (C ha ( ) A cr ) derived in Eqs. (23) and (25) with respect to the gridded dataset from Monfreda et al. 35 ); 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. 35 database and the overall methodology applied in the reconstruction of cropland/agricultural areas (as described above). 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 ha ( )

crops FAO
) and crop area provided in this study, as given in Eq. (38): The calculated ratio was then applied to the estimates of gridded annual crop specific harvested area (A crops (ha)) to obtain the crop harvested area harmonized at country level with FAOSTAT (C ha ( )

A crops
), 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):  www.nature.com/scientificdata www.nature.com/scientificdata/ (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 ha ( ) A cr ) 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).  Einarsson et al. 21 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. 35 , 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 1,30 . In particular, a recent study suggests that a large share of fertilizers are applied to permanent pastures in many European countries 21 . However, differences in the fractions of fertilizer applied to croplands and pastures have been reported as being a major source of uncertainty 5 . 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 41 . 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 fer Missing East EU (kg 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):  42 that provides global estimates of fertilizer production during 1925-1960 based on different data sources 43,44 . 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. 42  www.nature.com/scientificdata www.nature.com/scientificdata/ 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 crops Rate (kg ha −1 ) of crop harvested areas yr −1 ) and grassland (N fer grass Rate (kg ha −1 ) of grassland areas yr −1 ), as summarized in Eqs. (44)(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 crops Rate ) and grassland (N fer grass Rate ) 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)(47)(48)(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 (Nfer 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 crops Rate (kg ha −1 ) of crop harvested areas yr −1 ) and N fer grass Rate (kg ha −1 ) of grassland areas yr −1 ), respectively) to derive adjusted fertilizer application rates for crops and grassland (C N fer crops Rate (kg ha −1 ) of crop harvested areas yr −1 ) and C N fer grass Rate (kg ha −1 ) of grassland areas yr −1 ), respectively), as given from Eqs.   www.nature.com/scientificdata www.nature.com/scientificdata/ 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. 21 for the period 1961-2019. For the countries with incomplete record in Einarsson et al. 21 , 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 cr (kg yr −1 ) and N fer past (kg yr −1 ), 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:  For the second approach, we replaced the fertilizer application rates that are used in Eqs. (52) and (53) (C N fer crops Rate and C N fer grass Rate ) by the adjusted rates based on information provided by Einarsson et al. 21 . 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 5,27,28 .
To account for these uncertainties, we utilized information from two datasets of animal manure (i.e., FAOSTAT 45 and Einarsson et al. 21 ) 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) 46 . Einarsson et al. 21 , on the other hand, follows the approach of Lassaletta et al. 47 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. 21 , 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 app soil (kg yr −1 )) and manure left on pastures (Man left past (kg yr −1 )) for the period 1961-2019 from FAOSTAT 48 and Einarsson et al. 21 . 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. 29 that was derived using the spatial distribution of livestock counts from the Global Livestock Impact Mapping System (GLIMS) 49 and N excretion coefficients from the IPCC at a 5 arcmin spatial resolution for the time period 1860-2014. Zhang et al. 29 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 prod (kg ha −1 ) of grid area yr −1 )) of Zhang et al. 29 and the grid cell area (A grid (ha)) to calculate the country-level ratio of manure applied to soil to manure production (R app prod / soil (-)) and manure left on pastures to manure production (R left prod / past (-)), given below in the Eqs. (54)(55). Note that these ratios were estimated for each of two manure databases 21 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. 29 (Man prod (kg ha −1 ) of grid area yr −1 )), as presented in Eqs. (56)(57). Finally, we checked the consistency of gridded manure application rates to fall within a maximum permissible limit following Bouwman et al. 50 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. 21 data) between cropland and pasture, we used two different methodologies to account for uncertainties. First, following a similar approach as Tian et al. 51 , 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 cropland area (C A cr (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). 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. 21 data) according to the country-level proportion of manure applied to cropland and pasture derived from Einarsson et al. 21 . The latter study 21 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. 21

estimated ratios.
Overall, by utilizing two different input data sources of manure applied to soil and left on pasture (FAOSTAT and Einarsson et al. 21 ) 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 52 . 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) 53 . 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 (2022) 9:612 | https://doi.org/10.1038/s41597-022-01693-9 www.nature.com/scientificdata www.nature.com/scientificdata/ 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 54,55 . In this study, N fixation over cropland, pasture, forest and semi-natural vegetation (BNF cr , BNF past , BNF For and BNF NatVeg , 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)(64)(65)(66).   N fixation rates are derived for different crops and land use categories from previous studies 8,21,54 . For fodder crops, we calculated the N fixation rates using the fodder crop production and fixation dataset of Einarsson et al. 21 available at country-level during 1961-2019. The N fixation for the fodder crops in Einarsson et al. 21 was estimated assuming a linear relationship between the crop yield and N fixed by the crops following the approach of Lassaletta et al. 47 . 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). N fixation rates for pastures was taken as 5 kg of N ha −1 as suggested by Bouwman et al. 8 . 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. 54 . For semi-natural vegetation, we applied the same N fixation rate given for the natural grassland following Cleveland et al. 54 (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. Crop production. We obtained country-level datasets crop production from FAOSTAT for the period 1961-2019 56 . FAOSTAT provides production dataset for 173 crops and excludes dataset for fodder crops. For fodder crops, we derived the dataset from Einarsson et al. 21 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. 21 , 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 Wanmali www.nature.com/scientificdata www.nature.com/scientificdata/ (1984) 57 as cited in OWD 58 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 wheat,OWD (kg ha −1 yr −1 )) 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 wheat (kg yr −1 )) which is the product of the wheat yield (Yield wheat,OWD (kg ha −1 yr −1 )) and wheat harvested area (C A crops (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 wheat,OWD ) for time period 1850-1960 and the wheat production from FAOSTAT (Pro wheat FAO ) in the year 1961 by the crop production (Pro) of all other crops derived from FAOSTAT (for non-fodder crops) and Einarsson et al. 21  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 past (kg ha −1 yr −1 )) as harvested and grazed grass following a method proposed in previous studies 8,59 . 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 Rem past (-) with the N inputs to pastures (Inp past (kg ha −1 yr −1 )) adjusted for the losses of N after manure application to the soil (N losses (%)) as given in Eq. (72): The N removal coefficient and N losses for the animal manure are varying among studies 8,59 , 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 Rem past based on Bouwman et al. 8 and a value of 0.2 for N losses based on FAOSTAT. In the second scenario, we assumed a C Rem past 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. 60 and 0 for N losses .
N removal from forest. N removal from forest (Rem for (kg ha −1 yr −1 )) was estimated by adapting the method from Chang et al. 61 , 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 for (kg ha −1 yr −1 )) by the rate of forest N removal (C Rem For (-)) as given in Eq. (73).  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 (2022)  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 −1 yr −1 , 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 −1 yr −1 in the Brittany region in north-west France and 69 kg ha −1 yr −1 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 12,62 . Between 1930 and1960, 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 −1 yr −1 ), further increased by a factor of 1.9 between 1960-1990 (49 kg ha −1 yr −1 in 1990), and reduced by around 39% during 1990-2015 (30 kg ha −1 yr −1 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 −1 yr −1 ) for EU-28, Germany and the Danube river basin (Fig. 4a-c). Between 1950Between -1990, the values for the Danube river basin kept on increasing (up to almost 15 kg ha −1 yr −1 in the 1980s), EU-28 showed a more stable pattern (10-13 kg ha −1 yr −1 ). 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 −1 yr −1 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 −1 yr −1 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). www.nature.com/scientificdata www.nature.com/scientificdata/ 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 −1 of grid area yr −1 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 63 . 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.

Fig. 3
Total N surplus (kg ha −1 of physical area yr −1 ) 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.
www.nature.com/scientificdata www.nature.com/scientificdata/ • 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) 5 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 5 . To check for further consistency, we compared our N-budget with country specific estimates for Denmark 22 (1990-2010), the United Kingdom 23 (1990, 1995 and 2000-2018), France 24 (1940-2010), Poland 25 (1960-2000) and Germany 26,64 ,(1950-2017), while results for 25 European countries were additionally compared to estimates given by Leip et al. 20 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. 64 and Häussermann et al. 26 (for methodological details we refer to Ebeling et al. 65 ).   23 Poisvert et al. 24 N soil surface budget Agriculture 1940-2010 Country (FR) 24 Eriksson et al. 25 N soil surface budget Agriculture 1960-2000 Country (PO) 25 Häussermann et al. 26   www.nature.com/scientificdata www.nature.com/scientificdata/ 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 19 and Zhang et al. 5 . 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. 5  To determine the differences in the average N surplus from this study with respect to Eurostat data (respectively Zhang et al. 5 ), we calculated the country-wise relative difference Diff (%), as given in Eq. 74:  5 . 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 19 and for 70% of the countries with respect to Zhang et al. 5 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 19 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. 5 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. 5 , 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. 5 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 19 and Zhang et al. 5 . 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 19 and for 85% of the countries with respect to Zhang et al. 5 .
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. 20 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. 64 and Häussermann et al. 26 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 26,64 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 www.nature.com/scientificdata www.nature.com/scientificdata/ 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 66 relies on sales/trades data that mainly come from fertilizer companies of their member countries, FAOSTAT 41 collects information officially from its member countries via questionnaire and Eurostat 19 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 67 . 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. 5 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 28 . Another plausible cause of uncertainty in N budget is the different crop types considered in the studies, especially fodder crops 27 .
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 5,19,[21][22][23]25,26 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. 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).

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 5 , 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 40 ) using the spatial variability in land use areas (Ramankutty et al. 34 dataset), and crop specific harvested areas and crop production from the Monfreda et al. 35 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 www.nature.com/scientificdata www.nature.com/scientificdata/ 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_surplus.