Global impacts of nitrogen and phosphorus fertiliser use for major crops on aquatic biodiversity

The intensive application of nitrogen and phosphorus fertilisers on agricultural land to fertilise crops has caused eutrophication, the nutrient enrichment of waterbodies leading to excessive growth of algae, deoxygenation and loss of aquatic biodiversity. Life cycle impact assessments (LCIA) are often used to determine the environmental impacts of fertiliser use. However, the lack of suitable methodologies to estimate the fate and transport of nutrients from soils makes crop and regional impact comparisons challenging. Using a newly devised, spatially explicit nutrient fate and transport model (fate factor, FF) within an LCIA framework, this study estimates the global spatial-variability of nutrient loss from fertilisation of crops and their relative impact on aquatic biodiversity, specifically species richness. The newly devised FFs are based on the global spatially explicit nutrient model IMGE-GNM. The FF’s enable us to assess N and P’s fate and transport from indirect soil emissions (arable land, grassland and natural land) to freshwater environments. Additionally we improve the spatial resolution of existing soil FFs for N within marine environments from basin scale to 5 arcmin resolution. We applied our FF’s within current LCIA methodologies to assess the nutrient loading (midpoint indicator) and final aquatic biodiversity impact (endpoint indicator) from 17 crops. Our results identify strong variability in inputs, loadings and impacts due to differences in the fate, transport and impact of nutrients within the local environmental context. Such variability is translated into large differences between the popularly used nutrient use efficiency (NUE) indicator and final aquatic impacts caused by specific crops. Heavily produced crops (maize, rice, wheat, sugarcane and soybean) with the highest loading rates to receptors did not necessarily have the highest aquatic impacts. We identified rank variability exists at different metric stages (fertiliser inputs, receptor loadings, aquatic impacts) specifically for wheat and sugarcane. Our results showed high global spatial variability in aquatic biodiversity impacts with significant biodiversity loss outside of the highest production regions. Our study identified, global hotspots for biodiversity impacts depend on the local context that exist beyond the field (e.g. the fate and transport of nutrients to receptor environments, and the receptor environment's vulnerability). Aquatic impacts from fertiliser use for specific crop commodities should be considered in decision-making for strategic fertiliser pollution control and environmentally sustainable crop-commodity trade sourcing. The development of the improved FFs should be used to aid spatially explicit and site-specific LCIA nutrient studies from soils.


Introduction
Crop production to support a rapidly growing human population and demand per capita has led to intensive fertiliser use (Tilman et al. 2011). Fertilisers, namely nitrogen (N) and phosphorus (P), cause 78% of the global marine and freshwater eutrophication (Poore and Nemecek 2018), the nutrient enrichment of waterbodies causing excessive growth of algae, deoxygenation and biodiversity loss. Both fertiliser use and areas of eutrophication have grown rapidly since the 1950s and 1960s, respectively (Diaz and Rosenberg 2008;Morari et al. 2011). Identifying eutrophication hotspots and their sources from unsustainable fertiliser use, for specific crop commodities, is imperative to inform effective policies for reducing aquatic pollution from fertilisers and enable sustainable crop commodity trade financing and sourcing decisions.
Due to their extensive eutrophication impact, fertilisers have been typically modelled within life cycle impact assessments (LCIA) (Khan and Mohammad 2014). LCIA is one of the four main phases to a life cycle assessment (LCA), a popular tool for evaluating the environmental consequences caused by a product's life cycle (from resource extraction, through all subsequent processes involved in production, up to the end of the product's life). As a regional phenomenon, eutrophication is heavily dependent on the nutrients' transport, fate, exposure and effect in the local environment, rather than solely the fertiliser input and emission at the field (Rosenbaum 2018). The LCIA approach is to devise characterisation factors (CF) that may depend on the: • Nutrient proportion transferred through the environment and the residence time within a receptor (Azevedo et al. 2013a) (fate factor, FF) measured in units of time. The FF allows the calculation of a midpoint impact indicator representing a nutrient increase in water. • Exposure of the local environment to nutrient enrichment (exposure factor, XF) • Effect of nutrient loading on aquatic species (effect factor, EF). The EF allows for the calculation of an endpoint indicator relating to environmental damage.
The CF converts emission inventory data (collated at the inventory analysis phase of an LCA) into potential impact (I) by applying nutrient-specific (N or P) CFs for each elementary flow: Current LCIA methods differ in their emission inventory requirements, impact stage (midpoint or endpoint), geographical coverage, spatial resolution and emission pathways modelled. Hence, existing eutrophication cause and effect chains have variable levels of detail which have evolved over time. Previously, site generic indicators of freshwater and marine eutrophication typically occur at the midpoint between emission and environmental damage (endpoint). Midpoint indicators are used within several LCIA methods, such as ReCiPe 2008 (Goedkoop et al. 2008), EDIP 2003 (Hauschild and Potting 2004), IMPACT 2002+ (Jolliet et al. 2003) and CML 2002(Guinée and Lindeijer 2002. More recently, endpoint indicators have been incorporated into spatially explicit LCIA methodologies (ReCiPe 2016 (Huijbregts et al. 2016), LC-IMPACT (2016) and Impact World+ (Bulle et al. 2019)) that represent eutrophication damage to ecosystems. Even with current spatially explicit LCIA models identifying significant variability in the CF, most crop LCA databases and practitioners still use site generic models to estimate emissions of N and P fertilisers to waterbodies (Henryson et al. 2020;Notarnicola et al. 2017;Schmidt Rivera et al. 2017). Whilst site generic models simplify and reduce the data and modelling requirements for LCA practitioners, spatially explicit LCIA models offer opportunities to improve the quality and relevance of LCA studies, particularly with the growing demand to quantify environmental impacts from commodities to support sustainable trade and supply chain investing (Baitz et al. 2013;Patouillard et al. 2018).
Recent research on the CF has focused on the spatially explicit development of the FF, XF and EF (Azevedo et al. 2013b;Cosme et al. 2015;Cosme and Hauschild 2016;Payen et al. 2021). However, critical reviews of FF models have identified gaps in availability and spatial resolution of FFs for indirect emissions (e.g. fertiliser application to soils; non-point source) to surface waterbodies (Morelli et al. 2018). For indirect emissions, regionalised LCIA methods handle the transport of nutrients inconsistently and have inherent limitations (Morelli et al. 2018) (S.I 1 Part 1.1). These limitations include inconsistent methodologies (Helmes et al. 2012;Cosme 2018), coarse spatial resolutions (Payen et al. 2021;Helmes et al. 2012;Cosme 2018), FFs being estimated from globally unvalidated nutrient models (i.e. models not assessed against observed measurements) (Bulle et al. 2019) and lack of consideration for non-point source emissions to rivers (e.g. fertiliser application to soils) (Bulle et al. 2019). For marine environments, ReCiPe 2016 (Huijbregts et al. 2016) and recently LC IMPACT (2016) LCIA methodologies provide basin scale soil FFs for N, based on the Global NEWS2 Nutrient Model (Mayorga et al. 2010) through Cosme and Hauschild (2017). Current global freshwater LCIA FF methods developed by Helmes et al. (2012) are not based on validated nutrient models and only provide estimates from direct emissions to rivers (Helmes et al. 2012); they do not consider indirect soil emissions. Based on Helmes et al. (2012) methodology, ReCiPe 2016 assumes that 10% of P released to soil reaches freshwater, whilst IMPACT World+ (Bulle et al. 2019) makes no adjustments to FFs for indirect emissions. LC IMPACT (2016) updated FFs for P soil emissions using estimates from the SALCA model, through works by Scherer and Pfister (2015). Although the SALCA model was adjusted to the global scale (Scherer and Pfister 2015), it remains less accurate for regions outside of Switzerland, for which the model is programmed and validated to (Willmann et al. 2014). More recently, Payen et al. (2021) developed basin-scale FFs for N and P in freshwater using the Global NEWS 2 nutrient model. However, their results have not been implemented within any LCIA methodologies. Henderson et al. (2021) also developed fine resolution freshwater and marine eutrophication CFs using FFs by Helmes et al. (2012) and Cosme (2018) for the USA only. However, although the methodology is global in scope, only country level factors are provided at the global scale. Zhou et al. (2022) developed FFs for N in freshwater using IMAGE GNM (a global nutrient model) (Beusen et al. 2015) N emissions and loadings to freshwater as inputs. However, they did not provide FFs for P in freshwater or N in marine environments. Studies have identified that globally applicable and spatially explicit models for calculating soil emissions (i.e. leaching, erosion and runoff), which are more dependent on climate and soil conditions (e.g. temperature, pH, clay content and slope), are lacking in both N and P FF models (Notarnicola et al. 2017).
Few global LCIA studies exist for fertiliser use and their eutrophication impact. Scherer and Pfister (2015) showed large impact variability between crops for P use, with soybean causing the largest global impact to freshwater ecosystems per kilo produced. They additionally found that the highest emissions zones did not necessarily have the highest impacts on species richness. To our knowledge, global LCIA studies for N are lacking. However, nutrient transport and biophysical crop models have found high N and P emission hotspots in China, India, eastern USA, and central Europe (Liu et al. 2016). According to Diaz and Rosenburg (2008), these regions also show considerable marine eutrophication impacts. The nutrient use efficiency assessment has also been widely used to inform policymakers of regions where potential excessive fertiliser application is found (Brentrup and Lammel 2016;Zhang et al. 2015Zhang et al. , 2020. However, its use is limited, as it does not consider the environmental differences between regions or quantify nutrient losses to vulnerable receptors. Such limitations for LCIA fertiliser studies mean they have rarely been studied beyond the field edge and lack consideration for nutrient transportation to surface water environments and their associated eutrophication impacts (Schmidt Rivera et al. 2017;Liu et al. 2016Liu et al. , 2018. To summarise, current global LCIA studies, quantifying eutrophication impacts for crop products, have been limited by nutrient-type (Scherer and Pfister 2015), croptype (Cosme and Niero 2017) and/or spatial-resolution. To our knowledge, no study has utilised current spatially explicit LCIA methods to analyse the effects of N and P emissions from crops worldwide. In addition, current LCIA methodological limitations make it difficult to assess the emission, transport and fate of fertilisers from soils through the local environment. Hence this study's key contributions to the existing literature aims to provide the following: 1. To assess the impacts of eutrophication caused by N and P fertiliser use from crop commodities worldwide. We assessed impacts beyond the field edge to understand the aquatic biodiversity impact within receptor environments, which has rarely been studied in previous LCIA research. 2. To achieve objective 1 and fulfil current LCIA FF limitations, we developed a new spatially explicit FF model. This model fills gaps in current FF models (Morelli et al. 2018;Helmes et al. 2012;Cosme and Hauschild 2017), through developing spatially explicit, indirect soil emission FFs for N and P to fresh and marine water environments. Our FF model is based on concepts from a global nutrient model (IMAGE GNM) which has been validated against observed measurements. Additionally, our methodological framework is consistent across N and P to provide uniformity across FF methodologies, as recommended by many international initiatives, such as LEAP (FAO 2017) and UNEP/SETAC (Frischknecht et al. 2016).

Methodology
This study applies the four main phases of a standard LCA study to assess the eutrophication impact caused by 17 different crop products globally: • Phase 1: Goal and scoping-define the purpose and methodological specifications of the study. • Phase II: Inventory analysis-data collection and modelling of all inputs and outputs of the process making up a product's system. • Phase III: Impact assessment-translate inventory data into their contribution to the impact category. Quantify resource intensities and emissions that potentially cause environmental effects such as climate change, ozone depletion, acidification and eutrophication. • Phase IV: Interpretation-evaluating data quality and results to draw conclusions relating to the study's goal and scope.

Goal
The goal defined for this study was to determine the critical regional hotspots for eutrophication impacts caused by N and P fertiliser-use from 17 different crop commodities globally. The study evaluates fertilisers beyond the field edge, taking into consideration the importance of the local context (factors) that impact the removal of fertilisers and vulnerability of the environments they pass through. As such, we also compare our impact results to nutrient use efficient assessments conducted at the field. Due to data limitations and extensive requirements on data, we considered fertiliser emissions and aquatic biodiversity impacts associated with crop production for the farm stage only (i.e. we exclude other stages of food production along the supply chain that may be considered in studies focusing on a smaller area or particular commodity).

Eutrophication methodological scoping
The eutrophication methodological scoping for this study was done through a literature review. We wanted to identify the regional link between fertiliser emissions found at the field and the biodiversity loss in surface waterbodies through eutrophication. Current global LCIA methodologies still assume freshwater algae growth is limited by P and marine algae growth is limited by N (the limiting nutrient concept) (Finnveden and Potting 1999). Some recent studies have considered co-limiting N and P surface waterbodies; however, these are applied at national or regional scales (Henryson et al. 2018;Gallego Schmid et al. 2010). We adopted current global LCIA practices by applying the limiting nutrient concept due to the extensive literature review requirements and lack of adaptation of current global LCIA methodologies to cater for co-limiting or varying N and P-limiting surface waterbody assessment (Huijbregts et al. 2016;LC IMPACT 2016). For freshwater environments, we adopt the LC IMPACT (https:// lcimpact. eu/) methodology and for marine environments we adopt the methodology by Cosme and Hauschild (2017).
These LCIA approaches assess the nutrient loading (midpoint indicator) and final aquatic biodiversity impact (endpoint indicator). Nutrient loading represents the amount of nutrient reaching receptor environments and their duration of persistence, yielding a unit of kg ·yr. The final aquatic biodiversity impact is measured by the potential disappeared fraction of species (PDF). LC IMPACT quantifies the relative global species loss by putting regional species loss in perspective of the global species pool, leading to a globally normalised PDF of species (Verones et al. 2020). For freshwater, we solely consider fish species, as per LC IMPACT, whilst for marine environments, we consider demersal marine species (including fish, crustaceans, molluscs, echinoderms, annelids, and cnidarians), as per Cosme and Hauschild (2016). This study assumes global species pools for freshwater environments account for 15,000 fish species (LC IMPACT 2016), whilst marine environments account for 626 demersal (benthic and benthopelagic) fish and invertebrates species ). The adopted spatially explicit LCIA approaches use devised CFs (formed from the FF, XF and EF) to convert inventory fertiliser emissions into potential biodiversity impacts. For freshwater, impacts (I freshwater ) are estimated by: where E p reflects the P emissions (kg P ) and CF p (P characterisation factor) relates to the potential disappeared fraction (PDF) of species lost during a year due to 1kg of P loaded (PDF•yr•kg P −1 ). The CF p is defined by the newly devised FF p (fate factor, yr) and updated EF p (effect factor, PDF•kg P −1 ). FF p is the transported fraction of P from source to receptor and persistence in freshwater systems. EF p describes the potential decrease in fish species richness due to the P concentration in a specific freshwater type and ecoregion (Azevedo et al. 2013b). EF p is updated with improved hydrological data from PCR-GLOBWB 2 (Sutanudjaja et al. 2018). Spatially explicit XFs have not been devised for freshwater environments; only marine environments have spatially explicit XFs relating to oxygen depletion (Morelli et al. 2018). Further methodological details on the EF p is provided within the S.I 1 Part 2.2.
For marine environments, regionalised marine ecosystem potential impacts (I marine ) due to N emissions are estimated through marine eutrophication CFs (Cosme and Hauschild 2017): The E N (Kg N ) is the emissions of N. The CF N (N characterisation factor) is the potentially affected fraction (PAF) of species during a year (PAF yr kg N −1 ) defined by the updated FF N (yr) and the existing XF N (exposure factor, kg O 2 kg N −1 ) and EF N (PAF kg O 2 −1 ) (Cosme and Hauschild 2017). XF N and EF N vary globally between 66 large marine ecosystems (Cosme and Hauschild 2017). The improved FF N uses nutrient emissions fractions (f zN(i,j) ) from this study based from IMAGE-GNM (Beusen et al. 2015) with marine persistence rates taken from Cosme (2018). XF N reflects N's conversion to organic matter (phytoplankton biomass) to dissolved oxygen consumption in the bottom waters (Cosme et al. 2015). The EF N is the effect on biota due to dissolved oxygen depletion measured as the Potentially Affected Fraction (PAF) of species (Cosme and Hauschild 2016). To harmonise the endpoint scores in the LCIA framework, a conversion of the marine eutrophication endpoint CF from PAF to PDF is made using the conversion factor of 0.5 (Cosme and Hauschild 2016). The conversion factor assumes one-half of the species affected above their sensitivity to the hypoxia threshold (PAF metric) would disappear (PDF metric). This value is assumed based on the seasonality of planktonic production and biological processes, water temperature and stratification, and nutrient flows (Cosme and Hauschild 2016). The relative PDF metric was translated to an absolute metric through multiplying the PDF by species density distributions as per . We divided the demersal species densities by the global demersal species pool of 626 to obtain a global relative metric and harmonise our marine species impact results with our freshwater species result, following LC IMPACT (Verones et al. 2020).
We aggregated the biodiversity loss of freshwater and marine species to obtain the aquatic biodiversity loss found across receptors. To do so, we took the average PDF of species between marine and freshwater environments; this is one of the four aggregation methods suggested by Verones et al. (2015).

Fate factor methodological scoping
Our literature review showed that current FF models were not suitable for our spatial study regarding fertiliser impacts from specific crop products. Hence, we developed new FFs and applied them within our adopted LCIA methodologies. Here we develop new spatially explicit FFs that enable us to assess the fate and transport of N and P from indirect soil emissions (via surface runoff, erosion and leaching) from different landcovers (arable land, grassland CF N = FF N × XF N × EF N and natural land) to freshwater environments. Additionally, we improve the spatial resolution of existing soil FFs for N within marine environments (Cosme 2018) (Gleeson et al. 2014), Appendices 1 and 2). The updated datasets have a finer spatial resolution and reduce the number of assumptions within IMAGE-GNM, e.g. aquifer thickness. Each process in the development of the FF is driven by distinct local characteristics (e.g. slope, landcover, texture, temperature, soil loss, precipitation, soil drainage, soil organic carbon content) that alter the fate and transportation of nutrients. These processes have been described in detail by Beusen et al. (2015) for IMAGE-GNM.

Analysis scoping
Log regression analysis was used to identify whether the E, FF, XF or EF caused greater variability to the overall environmental impact found in fresh and marine waters as per Cosme and Hauschild (2017). To coincide with the previous LCIA literature (Payen et al. 2021;Cosme and Hauschild 2017), we assessed sensitivity through sensitivity ratios on primary, secondary and tertiary parameters that influence the final impact in marine and freshwater environments and build the FF, EF and XF as described by Strandesen et al. (2007). The sensitivity ratio is calculated as the ratio between the relative change in the model output and the relative change in the input parameters. We further assess some tertiary parameters to the FF which include arable land cover, temperature, precipitation, residence time, discharge and nitrogen load in rivers (Cosme and Hauschild 2017). These input variables were chosen based on their ability to be segregated easily within the model without impacting other input variables heavily. Annual temperature, precipitation and discharge between 2000 and 2015 were obtained from PCR-GLOBWB 2. We assessed the secondary parameter existing P soil content's impact on the FF using historical annual phosphorus soil contents between 1970 and 2000 which were obtained from IMAGE-GNM. For temperature, precipitation, discharge and P soil content, we varied each grid cell value by one standard deviation of the variable's annual distribution over the given time period. For all other variables, the input was varied by 10%.

Emission inventory analysis
Nutrient input and crop yield data were taken from http:// www. earth stat. org, as it provides data for 17 different crop commodities at the 5-arcmin resolution circa 2000 (Mueller et al. 2012;West et al. 2014;Monfreda et al. 2000). Although Earthstat provides the highest resolution available for global crop specific datasets, it is limited by data availability to 17 crops and to the year circa 2000 (1997-2003 period) solely.
Based on nutrient input data (Mueller et al. 2012;West et al. 2014) (kg N or P) for each crop type c, the soil nutrient balances for each 5-arcmin grid cell (i) is the difference between inputs and outputs of N and P: N inputs include biological N fixation (N fix ), atmospheric N deposition (N dep ), synthetic N fertiliser (N fert ) and animal manure (N man ). Outputs include N removed via crop uptake (N crop ) and N lost via volatilisation (N vol ). P has no fluxes with the atmosphere; therefore, only synthetic fertiliser (P fert ), manure (P man ) and phosphorus uptake (P crp ) are included. Positive soil balances give potential for nutrient to escape into the environment. Further information on the individual elements of the mass balance and their calculation can be found in the S.I 1 Part 2.1. Information on linking the emission inventory data to the FF can be found in the S.I 2 Part 2.

Results and discussion
Here our results and discussion present phases III and IV of an LCA study. Our results show: 1. The total aquatic impacts, nutrient loadings and inputs due to crop fertilisation. Unlike previous LCA studies, we identify eutrophication hotspots caused by the total production of crop commodities (i.e. not per unit mass or area but rather the total within a grid cell or region). Through this, we identify whether variability existed due to the local context or whether higher production and fertiliser inputs necessarily lend themselves to higher nutrient loadings and aquatic impacts. 2. The intensities of aquatic impacts, nutrient loadings and inputs per dry mass unit (tonne) of crop produced. This helps inform the trade-off between the amount of food produced and biodiversity impact, whilst fairly comparing performance across various crop types and regions. Although this is similar to previous LCA studies, here we assessed intensity by dry matter due to crops with high water contents, such as sugarbeet and potato, being disproportionately heavy. Assessment of intensities per dry mass makes crops more relatable to their end use as food, feed, fibre or fuel and support trade and supply chain risk analysis. 3. An assessment of our developed FF model to evaluate dominant pathways (runoff, erosion or leaching) for aquatic pollution/loading of nutrients, the influence of the local environmental context, sensitivity of distinct local characteristics and potential nutrient management strategies.

Crop specific
Maize, rice, wheat, sugarcane and soybean production account for 76% of the global dry mass for the 17 crops considered in this study (Table 1). Accordingly, fertiliser inputs (anthropogenic only) and nutrient loadings to aquatic receptors (freshwater and marine) for four of these five crops (excluding sugarcane) consistently ranked the highest against other crops. However, variability existed within (e.g. rice and wheat) and outside (e.g. potato) of these top five crops for their input, loading and impact rankings. Nutrient load to receptors and aquatic impact ranking showed a greater reliance on the fertiliser input than the total production (e.g. cotton and rapeseed). This suggests differences in crop fertilisation application and requirements influence the environmental impact more so than total production. Sugarcane, although ranked within the top five for production, ranked within the middle of the total crop distributions for fertiliser inputs, loads and impacts. For the top five impacting crops, sugarcane was replaced by barley which, whilst having 88% greater aquatic species impacts than sugarcane, only had 46% of Table 1 Global totals for production and the various metric stages (fertiliser inputs to fields, nutrient loadings to receptors and aquatic species impacts for both N and P). Global nutrient use efficiency (NUE) for each crop and ratios of load (R Load ) or impact (R Impact ) per kilo of fertiliser applied. Crops are ordered based on their total global production with ranking out of 17 ( sugarcane's dry mass but 35% greater fertiliser inputs. Sugarcane's total aquatic impacts were in a similar range to rye production, even though rye production is fifteen times less than sugarcane. Some crops showed large rank variability between different metric stages (fertiliser inputs, receptor loadings, aquatic impacts). Rye and sugarbeet had the highest raise in rank (6 rank difference) between nutrient input and aquatic species impacts, suggesting high nutrient delivery and vulnerability in their growing regions. Other crops showed a decline in rank from fertiliser inputs to aquatic impacts (e.g. cotton and groundnut) suggesting reduced nutrient delivery and low vulnerability within their growing regions. The least impactful crops, due to total production, included groundnut, cassava and oilpalm. Oilpalm's production had the least impact within the whole crop distribution and ranked 7th highest for dry production. We applied the nutrient use efficiency (NUE) indicator which is the ratio of nutrient applied to the crop uptake of nutrients. It is a well documented indicator of environmental sustainability (Brentrup and Lammel 2016;Zhang et al. 2015). Here the NUE is calculated using both fertiliser and natural nutrient inputs to understand the efficiency of fertiliser applied. Very low NUEs below 50% develop surplus N and P on fields, increasing the risk of nutrient delivery to vulnerable environments and economic costs (Brentrup and Lammel 2016). Thirteen of the 17 crops showed NUEs below 50% and included the top six most produced crops. Potato, sugarbeet and cassava had the worst NUEs below 20%. Optimal NUE's operating space is around 70-80%. Global groundnut and sunflower production were closest to this optimal range. NUEs above 90% (oilpalm) suggest soil mining and degradation potential due to insufficient nutrient inputs for crop production.
We compared crop ranking of the global NUE indicator to aquatic impact. The R ratio, which is the ratio between the fertiliser aquatic impact (R Impact ) or load (R Load ) to the fertiliser application, identifies the highest impact or load per unit of fertiliser applied to fields. Sugarbeet, rye, potato and soybean had the highest R Impact values suggesting the greatest inefficiencies, nutrient loadings and impacts to receptors are found for production of these crop commodities. Large ranking differences can be found between the NUE indicator and the R Impact for sugarcane, soybean, sorghum, millet and rye. Soybean and rye production are particularly impactful per unit of nutrient applied due to rye's vulnerability of receptors within its growing regions (seen by the R Impact value), soybean's high transport of nutrients to receptors in growing regions (seen by the R Load value) and soybean's high natural inputs through biofixation.
Our results identify strong variability between different metric stages (fertiliser inputs, receptor loadings, aquatic impacts) even when considering total production. This variability is translated into large differences between the NUE indicator, R Load and R Impact ratios at the global scale. The cause of this variability is due to differences in nutrient efficiencies on the field and variance in the fate, transport and impact of nutrients within the local environmental context (modelled through the various CF factors; S.I. 1 Part 3.3).

Eutrophication regional hotspots
Half of the world's total crop production considered in this study is found in North America, Europe, Central-Eastern South America, Eastern and Southern Asia. These regions have inherently much higher total fertiliser inputs than others, lending to substantial nutrient receptor loadings and considerable species richness impacts (Fig. 1a). The Yangtze, Mississippi and Ganges-Brahmaputra are specific basins within these regions with high biodiversity impacts from both N and P fertiliser use ( Fig. 1d and e). In general, biodiversity impacts were higher in marine environments than in freshwater environments and did not necessarily coincide spatially ( Fig. 1b-e).
Differences in the LCIA characterisation factor (CF) caused spatial variability in freshwater and marine environments (S.I Part 1 3.3 Fig 3). In marine environments, spatial variability was mainly caused by the farming processes and transportation of N from field to marine receptors (input-loading correlation; R 2 = 0.516; Fig. 2d), rather than the vulnerability of marine species (loading-impact correlation; R 2 = 0.99; Fig. 2d). Stronger correlations between P input-loading were observed in freshwater (R 2 = 0.754; Fig. 2e), however weaker loading-impact correlations (R 2 = 0.564; Fig. 1e) reflect that freshwater impact variability was influenced by fish species vulnerability. Hence, high freshwater species impacts were additionally found in basins with lower receptor loadings (e.g. Orange basin, Southern Fig. 1 Global maps of total aquatic biodiversity impacts from fertiliser use on all crops (a) assessed for nitrogen in marine environments (b) and phosphorus in freshwater environments (c). Each cell identifies the environmental impact to respective receptors (freshwater or marine) through fertiliser application at the cell's location. Total fertiliser inputs on the field (tonne), fertiliser loadings (tonne) to water bodies and impacts on aquatic species for nitrogen (d) and phosphorus (e) at the basin scale level are presented along with linear regressions and R 2 values. Names of basins refer to the top 5 basins with the most impact to species richness. Shape size represents the quantity of total crop production. Colours identify the biodiversity impact measured as the potentially disappeared fraction of species (PDF•yr); a commonly used LCIA biodiversity metric unit 18.2% and leaching dominated 67.4%. In P, surface runoff dominated 18.5%, erosion from new nutrient inputs dominated 0.5% and erosion from existing P in soils dominated 81.0% of the earth's arable land surfaces Africa). The Tigris-Euphrates basin is another clear example of this loading and impact variability, which identified as having the second highest loading rates for P, yet environmental impacts to fish species ranked 18th. Overall, low input-impact correlations in N (R 2 = 0.5; Fig. 1d) and P (R 2 = 0.607; Fig. 1e) meant high biodiversity impacts were additionally found outside of regions with the highest production and fertiliser application.

Local context and nutrient transportation; fate factor
The fate factor model includes three nutrient transport processes (surface runoff, erosion and leaching) that influence the emission of nutrients from fields to waterbodies. Figure 2 shows the spatially dominant emission pathways for N and P to freshwater bodies. Whilst both N and P are emitted via surface runoff and erosion, only N is considered to leach, as P is easily absorbed by soils (Beusen et al. 2015). As nutrients are applied to a field, surface runoff is considered to take a fraction of nutrient application (both natural and anthropogenic) before it has potential use for crop growth in IMAGE GNM. Surface runoff dominated 14.4% of the global arable land surface for N and 18.5% for P and was found mainly in North America, Central South America, Central Africa and South, East and South-East Asia.
Erosional processes for N and P are related to the quantity of soil loss per year in IMAGE GNM. N soil loss only takes account of existing N in the soil (i.e. without additional human or natural nutrient application) and is related to the amount of carbon found in soil loss. Therefore, it is unrelated to the quantity of fertiliser applied to fields. In P, erosional processes are related to the mass balance from new P inputs, historical P inputs and existing P found in the soil. Consequently, for both N and P, an inherent quantity of nutrient is lost via erosion from arable fields that is uninfluenced by the quantity of N and P fertilisers applied on the field. The inherent quantity of nutrient loss via erosion for N dominated 18.2% of the earth's surface, whilst for P, inherent soil loss dominated 81.0% of arable land surface and 0.5% was dominated by new P inputs to soil loss.
Existing erosion was the primary cause of P removal from agricultural soils and leaching was the primary cause of N removal. N leaching dominated 67.4% of arable grid cells and existed in most major agricultural basins (e.g. the Mississippi, Ganges-Brahmaputra, Danube and Yangtze). N fertiliser is leached into the groundwater where it can persist over many years before leaching into rivers. Although IMAGE GNM models nitrogen persistence in groundwater up to 50 years, this could be much longer due to storage within the vadose zone, which is not modelled within IMAGE GNM (Ascott et al. 2017). Ascott et al. (2017) showed that global storage volumes of nitrate in the vadose zone are between 605 and 1814 teragrams. Storage of nitrate in the vadose zone may delay the transient delivery of nitrates through groundwater to rivers, thereby delaying eutrophication impacts. The uncertainty in the quantity and length of storage time within the vadose zone and groundwater can create negative implications on nutrient management strategies applied to fields or receptors as leaching rates into river receptors are dependent on the historical fertiliser inputs, the quantity of denitrification in soils and riparian zone and the length of time of storage in the vadose and groundwater. Globally variable historical fertiliser use means that whilst some regions have increased their fertiliser use, others have decreased use over time. Here we factor in transient (variable) historical fertiliser use by developing a transient state (TS) versus steady state (SS, constant input of fertiliser over time) factor which increases, or decreases, N loads for the year 2000, depending on the historical load existing within aquifers. We applied a maximum TS/SS factor of 2 due to high uncertainties found within groundwater datasets (Sutanudjaja et al. 2018).
For more details on the FF model methodology, validation, results and comparisons with other FF models, we direct the reader to the supplementary material S.I 2: Development of a Spatially Explicit Fate Factor Model at the 5 Arcmin Resolution. The FF results for all land covers (arable, grassland and natural) are found in this supplementary material.

Crop specific
To fairly compare performance across various crop types, we assessed the inputs, loadings and impacts of crop production per dry mass unit (tonne) of crop produced (Fig. 3). The FF model suggests there is an inherent loss of nutrients from fields without additional anthropogenic fertiliser input. Hence, we separated natural arable inputs, loads and impacts from those associated with anthropogenic fertilisation inputs. Whilst natural inputs include N deposition and fixation inputs, natural arable loadings include loadings from natural inputs and N and P erosion from arable soils (using natural landcover as a baseline for erosion; see S.I Part 2 2.4). Separating fertilisation flows from natural flows allows us to assess the impact of anthropogenic fertiliser within any given year. Crops are ordered based on their fertiliser impact to fertiliser input ratio (R Impact , Table 1). Across aquatic environments (marine and freshwater), the potential disappeared fraction (PDF) of species was highest for rapeseed, potato and rye. Rapeseed impacts from fertiliser use were over four times higher than the mean impact from fertilisers across crop commodities (2.78 × 10 −12 PDF•yr/tonne crop). The lowest impacting crops included cassava, sugarcane and oilpalm and were related to low fertiliser inputs per tonne of dry crop. Fertilisation inputs and receptor fertiliser loadings for rapeseed were the highest across all crop distributions. Rapeseed fertilisation rates were sixteen times greater than the lowest fertilised crop, oilpalm. This resulted in loading rates over 40 times higher for rapeseed than oilpalm. Rapeseed is an oil crop; therefore, switching from rapeseed to another oil crop (e.g. soybean, sunflower or oilpalm) may reduce global eutrophication impacts. However, it is important to note that increasing the production of another crop may exacerbate their eutrophication impacts or have other environmental consequences such as deforestation and terrestrial biodiversity loss (Qaim et al. 2020).
Soybean's atmospheric nitrogen fixation abilities increased soybean's natural N inputs, which had over 10 times greater natural inputs and more than double the quantity of nutrient loadings (not considering fertiliser loadings) to receptors, than any other crop. However, fertiliser loading rates and overall marine impacts of N for soybean were low compared to N inputs, suggesting soybean benefited from strong nutrient use efficiencies in N found on the field, low N nutrient delivery to marine receptors and low marine species vulnerability in soybean's growing region. Opposingly, and coinciding with previous literature, soybean had the highest freshwater fish species impacts due to high vulnerability of freshwater receptors within soybean's growing regions (Scherer and Pfister 2015) (S.I. 1 Part 3.1). This indicates nutrient management strategies should be directed towards protecting freshwater environments for soybean production.
N fertiliser application rates were on average five times greater than P fertiliser application rates and N loading to marine environments was an order of magnitude or two greater than P loading to freshwater (S.I. 1 Part 3.1). Consequently, marine species impacts were greater than freshwater fish species impacts. For example, rye and rapeseed caused 104-and 116-times greater impacts in marine environments compared to freshwater environments, respectively. This was due to greater nutrient loading rates to marine environments Fig. 3 Global production weighted averages of N and P fertiliser field application (blue), aquatic receptor loadings (green) and aquatic impacts (red) per tonne of crop for seventeen different crops. Crops are ordered ascending based on their ratio of fertiliser aquatic impact to fertiliser inputs (R Impact ); high ratios mean better use of nutrients on the field for crop production and reduced nutrient transport to aquatic receptors, low ratios suggest inefficient use of nutrients on the field and increased nutrient transport to receptors. Separate graphs for N and P are found in the S.I. 1 Part 3.1 and higher marine vulnerability where rye and rapeseed are grown, specifically in Europe. Whilst rapeseed, potato and rye identified as having the highest impacts to marine species, soybean, rice and potato showed the greatest impacts to freshwater fish species. Additionally, greater variability existed in freshwater impacts between loadings and impacts (e.g. groundnuts versus barley), indicating strong spatial variations in freshwater fish species vulnerability in specific growing regions. For P fertilisation, we identified multiple cropping systems having low R Impact values given their high P fertilisation inputs (e.g. sorghum, sunflower, millet and rapeseed). This is due to lower nutrient delivery and, in some cases, lower freshwater receptor vulnerability, e.g. for rapeseed in Europe.

Country level performance
At the country level, we assessed the performance for specific crop commodities per unit of dry mass produced and applied the NUE indicator. We assessed the relationships between production, NUE, loading and impact metrics for all crops and countries. We show the relationships which were significant (P > 0.05) in Table 2.
One would expect countries with high production to benefit from economies of scale and have suitable NUEs that translate into lower impact intensities within waterbody receptors. However, at the country scale, we found no significant relationship existed between production quantities and intensities of inputs, loading or aquatic impacts. Similarly, no strong relationships were identified for the combined NUE indicator with loading or impact metrics. For N and P individually, no strong relationship (< −0.2 or > 0.2) existed between receptor impacts and the NUE indicator (−0.11 for N and 0.15 for P). The lack of a strong correlation between the NUE indicator and aquatic impact metrics makes the NUE indicator unsuitable to assess environmental sustainability from fertiliser use at the country scale. This solidifies the need for endpoint indicators that measure biodiversity impacts rather than indicators related to resource use efficiencies or intensities.
Strong positive relationships were identified between aquatic impacts and combined loadings (0.59). However, considering N and P individually, a notably stronger relationship was identified in marine environments (0.67) than in freshwater environments (0.19), indicating strong variability in freshwater vulnerability. For N and P in freshwater, a strong positive relationship between N and P loading existed (0.32). However, only N freshwater loading showed significant relationships with N marine loading (0.37). Our results show strong variability between all metrics (production, NUE, loading and impact metrics) suggesting the local Table 2 Correlation matrix of production, NUEs and intensities for combined N and P and separate at the country level on the global scale for relationships which were significant (P > 0.05).
Strong correlations are highlight bold (< −0.2 or > 0.2) Combined N and P Aquatic environments N to marine and freshwater environments P to freshwater environments For each crop, country level scatterplots using production, NUE, loading and impact metrics were devised. Although the results for all 17 crops can be found in the S.I 1 Part 3.2, we present here (Fig. 4) for illustrative purposes, the results for the two most heavily produced crops worldwide, maize and rice. The impact presented here is the overall impact including the underlying natural impacts measured in PDF•yr/tonne dry crop. We assessed the global median rather than mean due to the strong positive skewness of the data.

Production NUE N and P Load N and P Aquatic species impact NUE N N load marine N load freshwater N species richness NUE P P load
At the country level, the global combined median NUE, for both N and P, was below 50% (48.65%). This suggests surplus N and P on fields increases the risk of nutrient delivery to vulnerable environments and economic costs. Most crops showed country level NUEs above and below the global median. However, some specific crops showed almost all countries performed with NUEs below this value (e.g. potatoes) or above (e.g. oilpalm). Many African countries had combined NUEs well above 100% such as Madagascar (oilpalm NUE = 3278%; cotton 578%), Niger (rice NUE = 1829%) and Ghana (Groundnuts NUE = 289%). The NUE indicator is a suitable indicator to understand resource use and efficiencies found on the field. However, its use to assess Fig. 4 Intensities of nitrogen and phosphorus loadings and impacts to waterbody receptors per unit ton of yield at the country level. Countries with the smallest production (for a total of 0.01% of global total production of each crop) are discarded; shapes represent different continents; sizes represent production of each country; colour represents NUE (%); dashed red line represents the global mean of loadings and impacts emissions and environmental impact is limited as the NUE indicator does not consider the inherent loss of nutrients from soils, the reduction of nutrients during their delivery to receptor environments from fields and the varying level of biodiversity vulnerability found within receptor environments.
In marine environments, the greatest impact to marine demersal species, on orders of 3.00 × 10 −10 to 1.00 × 10 −11 PDF•yr/tonne, was mainly found in East Asia by multiple crops (e.g. soybean, rye, rice and cotton) and in Europe along the North Sea by rapeseed, potato and rye production in countries such as Latvia, Sweden and Belgium. Most NUEs for N in these regions were below 50%, suggesting poor nutrient use efficiencies are found for multiple cropping systems. However, additionally strong nutrient delivery and receptor vulnerability play major roles in increasing the environmental impact in these regions. One example is the Baltic Sea in Europe which has a residence time of 21.83 years, 9.4 times greater than the average residence time found globally in large marine ecosystems (Cosme and Hauschild 2017). The residence time increases the nutrient persistence in receptors, increasing eutrophication potential. Additionally, in East Asia the seas surrounding Korea have high marine vulnerability (S.I. 1 Part 3.3). For nitrogen loading to marine environments, the highest loadings were found in Europe (e.g. Bulgaria, Turkey, Italy and Latvia) and South Asia (Bhutan and Bangladesh) for mainly rapeseed and soybean production. The highest loading to marine environments was by cotton production in Bulgaria, where 41.93 kg N/tonne dry crop reaches marine environments causing high marine impacts (4.84 × 10 −11 PDF•yr/tonne). The lowest marine impacts were mainly found in Africa, in particular Uganda, Niger and Mali, on orders of 1.00 × 10 −15 to 2 × 10 −18 PDF•yr/ tonne, for multiple production systems such as millets, maize, sugarcane and cassava. This is mainly due to countries being inland where emissions are heavily reduced, through retention in soils and river systems, before reaching the coastal environments and low vulnerability of marine receptors surrounding Africa. Countries in Central Asia such as Turkmenistan, Azerbaijan and Uzbekistan showed zero impact to marine species, as hydrological flows belonging to these countries are not connected to large marine ecosystems.
Unlike marine impacts, freshwater fish species impacts in Northern Europe were some of the lowest, particularly in the Netherlands, Denmark, Sweden and Belgium, for crops such as barley, wheat, sugarbeet and potato. Low NUEs in P suggest the low freshwater fish species vulnerability (characterised through the EF) plays a major part in reducing the freshwater impact (S.I. 1 Part 3.3). The highest freshwater species impacts were found in Lebanon, Israel and Lesotho, on the order of 3.5 × 10 −10 to 1.03 × 10 −10 PDF•yr/tonne. These countries are all semi-arid, sub-humid climates and hence have increased receptor vulnerability. The vulnerability of freshwater fish species (characterised through the EF for P) is dependent on the fish species richness, climate, volume of water and waterbody type (river or lake/reservoir) found within ecoregions; see S.I. 1 Part 2.2 for the EF methodology and LC impact for further information. As such, ecoregions with a dry subhumid or semi-arid climate, with lower water volumes and greater quantities of lakes and reservoirs, increase fish species densities and their inherent vulnerability within water-scarce regions.
Our results suggest fertiliser inputs for specific crops vary across countries and production scales. The local context (fate and transport of nutrient and vulnerability of receptor environments) plays a great part in enhancing or diminishing environmental impacts from fertiliser use, suggesting a strong reliance on regionality of environmental impacts rather than quantity of fertiliser use, or type of crop grown.

Global regression analysis and sensitivity analysis
A spatial regression analysis was used to identify whether the emission inventory (E), fate and transport factor (FF), exposure factor (XF, in marine environments only) or species vulnerability in aquatic receptors (EF) caused greater variability to the overall impact (I) caused by the production of wheat globally (S.I 1 Part 3.5). We chose wheat as it is the most heavily produced crop globally with the greatest spatial coverage. The analysis showed that whilst the intrinsic spatial variability of the CF within marine environments was primarily driven by the FF, in freshwater environments the FF and EF had equal influence. The emission inventory had a low influence for both N and P. Drivers of spatial variability are indicated through high coefficients of determination (R 2 ), relatively low sum of squares, mean square and standard error. The LCIA impact study sensitivity to the four main input parameters, i.e. E, FF, XF and EF, was assessed through sensitivity ratios (SR) to coincide with previous LCIA literature (Cosme and Hauschild 2017) (S.I 1 Part 3.6). As the impact found in marine or freshwater environments is a linear multiplication function, each of the primary parameters showed an expected SR = 1. We assessed sensitivity to primary parameters that build the FF, EF and XF. The FF model showed neutral sensitivity ratios (SR = 1.00) to the soil emission fraction along all emission pathways. For the downstream retention reduction factor, strong positive SR ratios suggest increasing retention in streams increases the FF value by over 5 times in freshwater and 9 times in marine environments. The persistence rate showed an almost equal negative SR (SR = −1.1). The EF in P showed neutral sensitivity (SR = 1.00) to fish species richness within ecoregions. However, the volume of water (all rivers, lakes and reservoirs) within an ecoregion had an average negative SR of −0.91 (S.I 1 Part 3.6). This suggests freshwater-scarce regions such as Lebanon have increased fish species vulnerability. Lebanon had the greatest impact to freshwater fish species richness per tonne of groundnut, potato, sugarbeet and wheat produced globally (S.I 1 Part 3.2). Cosme et al. (2015) and Cosme and Hauschild (2016) studied the sensitivity of the XF and EF to parameters in marine environments respectively. They found LMEdependent primary production rates (photosynthesis rate to convert nutrients to organic plant mass), secondary producer assimilation fraction (transfer of organic plant material grazed by other trophic levels e.g. zooplankton) and primary producer sinking fractions (plant material sinking to the ocean floor) had the highest positive SRs contributing to the XF in marine environments (Cosme et al. 2015). For the EF, species sensitivity to hypoxia (as a lowest-observedeffect-concentration) showed little sensitivity (Cosme and Hauschild 2016).
We additionally assessed secondary and tertiary parameters that build the primary parameters for the FF model (S.I 1 Part 3.6). For secondary parameters, we assessed SRs for different emission pathway fractions. The greatest sensitivity was identified in the erodibility of P in soils (SR = 0.99) and for N the greatest sensitivity existed for leaching (SR = 0.37-0.54) and soil erodibility (SR = 57-0.72).
For tertiary parameters, we split parameters based on their impact to soil or water compartment of the FF. For the FF soil compartment, average SRs ranged between −0.03 and 0.15. Suggesting small changes to the FF occurred from changing individual parameters such as temperature, land cover and groundwater depth. Whilst temperature in the soil compartment showed a small positive ratio (0.016-0.13) for its impact to the leach fraction, temperature showed strong negative SRs in the water compartment of the FF, with an average range across receptors from −3.97 to −0.83. The highest sensitivity to temperature were found along a band in the Northern Hemisphere, crossing over Canada and Russia, where some of the highest FFs currently exist. Studies have shown Canada and Russia will be prime agricultural countries under the impacts of climate change in years to come (Frischknecht et al. 2016). Although a rise in temperature may cause fate and transportation of N and P to waterbodies to decrease (due to an increase of nutrient uptake in waterbodies and increased N losses to the atmosphere through denitrification in the soils, riparian zone and waterbodies), a rise in temperature may have an adverse effect on the ecosystem's vulnerability and other environmental impacts.
Under the influence of climate change, clear trade-offs exist between nutrient runoff reductions and increasing greenhouse gas emissions which develop positive climate change feedback loops. Opposingly, temperature and other climatic factors (rainfall) may also impact the vulnerability of the waterbody receptor, through adverse variations in exposure and effect factors, that may increase the overall impact on species richness. Discharge showed a strong positive SR impact in marine environments whilst having negative SRs in freshwater environments, identifying the strong influence of advection to increase the removal of nutrients from freshwater systems and loading to marine environments. Hence, obstructions to river flow paths, such as dams, significantly impact the FF in marine environments. The hydraulic load and residence time also had a strong influence in marine environments. Whilst increasing the hydraulic load reduced retention which increased nutrient delivery to marine environments, reducing the residence time reduced the nutrients persistence in freshwater environments thereby increasing nutrient delivery downstream. Overall, marine environments showed greater sensitivity to tertiary parameters. We urge caution around the sensitivity of results, as input parameters have strong influences on other parameters which would vary in parallel, particularly hydrological parameters.

Nutrient management strategies
Under the assumption of the limiting nutrient concept (Harris 1986), we assessed the inflection point used by Helmes et al. (2012) to detect the difference between longer and shorter FF persistence values for arable soils. For arable land covers, a persistence exceeding 0.08 and 0.11 days for P and N reflected longer persistence in freshwater environments respectively, whilst in marine environments persistence inflection increased to 5.06 days (S.I. 1 Part 3.7). Persistence hotspots for transportation of nutrients (from all land covers) were in every continent and depended on the geomorphology and hydrological characteristics. Using Europe as a case study example, we identified the spatial location of lower and higher persistence values based on the inflection point (S.I. 1 Part 3.8). Coastal regions showed lower persistence for freshwater environments and higher for marine environments as nutrients are quickly transported out of the river system. Freshwater FFs were higher inland, within the upper reaches of river networks, as nutrients travel through the river network over the longest distance. Mountainous regions such as in central Russia, west of Lake Baikal, showed particularly higher FF values (FFs > 100 days) (S.I 2 Part 2.2 Fig  5). However, emissions further from the river system flow paths have lower persistence values, caused by the removal processes (denitrification, nutrient uptake and advection) and reduced soil emissions (surface runoff, erosion, leaching), meaning less nutrients reach main river systems. River systems with dams (e.g. the Volta River in Ghana) increased freshwater FFs upstream of the dam and decreased FFs to marine environments due to small hydraulic loads (dependant on high residence times) passing through dams, increasing retention. In the extreme case of the Volta River in Ghana, the particularly small hydraulic load (Akosombo Dam = 5.44 m/yr < Global median = 153.7 m/yr) passing through the Akosombo Dam entirely retained N nutrient, inhibiting delivery to marine environments. Hence, agricultural activities upstream of lakes and reservoirs may cause higher eutrophication impacts to freshwater environments due to the longer persistence of nutrients within lakes and reservoirs. Downstream of lakes and reservoirs, nutrients may persist for short times in the freshwater systems. However, this may increase the nutrient delivery to marine environments where they may persist for longer periods, causing greater marine eutrophication impacts.
Lower inflections were identified when assessing FF persistence values from grassland and natural landcovers where arable land existed. FF differences between different land uses are constant at the grid cell level. For erosion, arable and grassland FFs are 45.3046 and 2.4138 times greater than natural land erosion for the same grid cell (Beusen et al. 2015). Similarly for surface runoff and leaching, arable and grassland FFs are 8 and 2 times greater than natural lands (Beusen et al. 2015). This suggests land conversion to natural or grass landcovers support reducing the overall fate and transportation of nutrients to receiving waterbodies. A comparison of our FFs presented here has been made with previous studies (S.I 2 Part 2.3). Overall, our results, as well as previous studies estimations, suggest assuming diffusive emissions are 10% of direct emission is not sufficiently robust enough to estimate diffusive emission FFs. Additionally, large FF variations exist between studies, identifying the strong reliance of FF models on robust global nutrient models.
Going beyond the limiting nutrient concept (i.e. ignoring the limiting nutrient concept), we assessed the potential algae production in freshwater from N and P using the Redfield Ratio (Hecky et al. 1993). Global total potential algae production from N in freshwater equated to 8.66 Gg Algae /yr, and for P 4.23 Gg Algae /yr, in marine environments potential algae production from N equated to 52.7Gg Algae /yr. Specific hotspots for N potential algae production were in Eastern North America, Europe, Eastern and South East Asia and Central South America (S.I 1 Part 3.9). This approach does not accurately represent global algae production as freshwater environments are generally P limited. However, research has shown freshwater and marine environments can display N limitation, P limitation or co-limitation and the limiting nutrient could change both seasonally and spatially (Ngatia et al. 2019;Paerl 2009). Hence, a spatially explicit approach should be used for catchments that are co-limiting such as Lake Taupo, New Zealand (Pearson et al. 2016). Nevertheless, the indication of large freshwater N loadings, causing detrimental water quality and potential algae growth, suggests the need for joint nutrient management strategies.

Comparison with other studies
Few global LCIA studies are comparable with our study. This is because our study is the first (to our knowledge) to use soil FFs for both N and P on individual crops. Thus, we estimated the emission quantity at the edge of the field to make our study more comparable with previous research.
For N, there are no comparable global scale LCIA studies. Huang et al. (2017) evaluated species richness impacts from crop production at China's provincial level, which were not directly comparable to this study. We therefore compared our study to large scale nitrogen loss models (Liu et al. 2010(Liu et al. , 2016Bouwman et al. 2009;2013;Mekonnen and Hoekstra 2015;Sutton et al. 2013). The ratio of total N inputs to total N outputs (0.33) at the field were within range of previous estimates (0.22-0.55) and were most similar to Beusen et al. (2015) (Table 3A). For total P input and output ratios, our result (0.16) is within range of previous studies (0.025-0.21) (Table 3B). A high input loss ratio is the result of accounting for the initial P content in soils, which considers historical fertiliser input, when losses are calculated. Additionally, our initial mass balance and fertiliser inputs are much smaller than previous studies. Furthermore, crop nutrient residuals were not considered within our mass balance unlike in other studies.
For P, our emission rates are three to ten times greater than previous estimates for specific crops (Table 3C) (Scherer and Pfister 2015) and an order of magnitude greater than the Ecoinvent datasets. The impacts to species richness in freshwater environments were 2 orders of magnitude greater in this study than Scherer and Pfister (2015). Scherer and Pfister (2015) identified soybean production as the single most impacting crop per mass of crop on relative species richness. Here we found soybean impacts to be mild compared to other crops using Scherer and Pfister (2015), using the EF factor taken from Azevedo et al. (2013b) which does not consider the spatial variability of fish species densities. In this study, LC IMPACT method incorporates fish species densities within the EF factor which increased vulnerability in soybean production regions; ranking soybean the worst impacting crop from P fertilisation.

Limitations and uncertainties
Here we use annual mass balance approaches to quantify the emission inventory. This is to harmonise LCIA methodologies with IMAGE-GNM modelling and current data limitations. Although mass balance approaches are popular methods (Bouwman et al. 2009;Sciences and of. 2013;Liu et al. 2010), biophysical crop-growth models (designed for assessing crop yields) have been used with assumptions for fertilisation schedules to quantify nutrient emissions at the field (Liu et al. 2016(Liu et al. , 2018. Liu et al. (2016) identified three commonly used fertilisation schedules predicting similar nutrient losses, with only minor differences in yield output, and recorded similar losses to annual mass balance approaches. Biophysical crop models may provide more robust on-field results but data limitations at a global Table 3 Comparison of estimated global total nitrogen (A) and phosphorus (B) emissions at the edge of the field, and ratio of total losses to total inputs with results from previous studies. Comparison of phosphorus inputs, emissions and impact intensities (C) for the crops maize (M), rice (R), soybean (S) and wheat (W) with previous literature Sutton et al. scale (including, nutrient type applied, application schedules, management practices and technology used) mean differences in nutrient losses are insignificant (Liu et al. 2016(Liu et al. , 2018. As IMAGE-GNM has conducted a full sensitivity analysis on the three main components of the model (nutrient delivery, retention and export), it is beyond the scope of this study to undertake a full parameter sensitivity and uncertainty analysis of the FF. Here we provided a sensitivity analysis on the primary, secondary and some tertiary factors for the FF. IMAGE-GNM was validated by observed concentration data found in the Mississippi, Rhine and Meuse Rivers, as well as a further 125 rivers in Europe, for which they identified model performance as acceptable. For some of the observed monitoring data stations, poor agreement was identified and associated with peak periods of high rainfall. Beusen et al. Temporal differences in emission pathways (e.g. surface water versus groundwater) are not defined within the FF; this is commonly not accounted for within LCIA approaches (Helmes et al. 2012;Cosme and Hauschild 2017;Scherer and Pfister 2015). Furthermore, the initial FFs are devised for the year 2000 and therefore may misrepresent any other given year. Improving the datasets used within IMAGE-GNM equations could provide annual FFs between 2000 and 2015, in line with the updated hydrology model PCR-GLOBWB 2 (Sutanudjaja et al. 2018). Annual changes to environmental variables such as temperature, precipitation, discharge and soil nutrient content would impact soil loss, leaching, denitrification rates, sedimentation, nutrient uptake and residence times, all of which impact the FF. This would provide significant improvement to the development of inter-annual variations of FFs within LCIA methodologies. Developing FFs at other annual timescales does not benefit this study where indicator datasets are circa 2000. Crop-specific gridded data on yield and harvested area, on which the fertiliser and emissions estimates depend on, is currently only publicly available for 1997-2003 period (circa 2000). Furthermore, current eutrophication LCIA methodologies relating to the XF and EF do not provide inter-annual temporal differences within their factors.
The modelling presented in this study and IMAGE GNM do not fully consider the positive or negative impacts of specific farming techniques that may reduce both inherent nutrient loss and additional fertiliser loss from arable soils. However, the inclusion of farming or remediation techniques along FF nutrient emission pathways could be formed at the grid cell level through applying reduction factors for effective agricultural nutrient management practices. To devise such factors, further research into suitable reduction factors for farming and remediation techniques are required, as well as further data collection on the effectiveness of nutrient remediation practices found globally. Farming practices such as terracing, contour cropping and no-till planting are all used to reduce soil loss and runoff but may increase leaching processes. Applying factors for these farming techniques could be applied to future LCA studies allowing the FF model to be more practical and encourage nutrient loss management to leverage FF results at the grid cell level.

Conclusion
Here we presented an LCA study to assess the eutrophication impact in freshwater and marine environments caused by specific crop fertiliser applications globally. To assess the fate and transport of fertilisers from field to aquatic receptors, we developed a new spatially-explicit (5 arcmin) FF that could be applied within existing LCIA frameworks. The steps to create a new FF were necessary due to the limitations existing FFs have for nutrient movement from agricultural soils.
Our study goes beyond the nutrient excess found at the field by assessing the nutrient transportation to aquatic receptors and quantifying aquatic biodiversity impacts. Overall, our results identify strong variability at different metric stages (fertiliser inputs, receptor loadings, aquatic impacts) even when considering total production. Although crops with the highest total production tended to have the highest aquatic impacts (e.g. maize, rice and wheat), we found differences in total crop fertilisation inputs influenced the environmental impact more so than total production (e.g. sugarcane versus barley). The variability at different metric stages (particularly the strong variance in the fate, transport and impact of nutrients within the local environmental context) goes beyond the popularly used nutrient use efficiency (NUE) indicator's capabilities. This suggests the NUE is not a suitable indicator to define environmental sustainability and impacts caused by fertiliser use for specific crop commodities. Whilst the NUE indicator alone identified potato, sugarbeet and cassava with the worst global NUEs below 20%, our study's ratio of aquatic impact to fertiliser input ratio (R Impact ) identified sugarbeet, rye, potato and soybean as having the greatest inefficiencies, nutrient loadings and impacts to receptors per kg of fertiliser applied.
Regions with intensive agricultural production (Yangtze, Mississippi and Ganges-Brahmaputra) inherently had much higher total fertiliser inputs, substantial nutrient loadings and considerable species richness impacts. However, outside of basins with the highest production, differences in the LCIA characterisation factor (CF) caused spatial variability in freshwater and marine environments. Whilst spatial variability in marine environments was primarily caused by the farming processes and N delivery to marine environments, spatial variations in freshwater fish species vulnerability had a greater influence. High marine impacts were identified in Europe, specifically where rye and rapeseed are grown. In freshwater, ecoregions with a dry subhumid or semi-arid climate, with lower water volumes and greater quantities of lakes and reservoirs, have higher fish species densities and inherently high vulnerability within water-scarce regions (e.g. Lebanon and Greece).
Considering the performance of crop production (i.e. intensities per unit of dry mass crop), rapeseed, potato and rye had the greatest aquatic impacts whilst cassava, sugarcane and oil palm had the least. Using rapeseed alternatives (soybean, sunflower or oil palm) may reduce global eutrophication impacts. Due to soybean's atmospheric nitrogen fixation abilities, soybean had the highest natural inputs and loadings, but low marine species impacts due to low marine vulnerability within growing regions. However, high vulnerability of freshwater fish species in soybean growing regions meant soybean had the highest impacts from P fertilisation across all crops, indicating nutrient management strategies for soybean should be directed towards protecting freshwater environments. N fertiliser application rates were on average five times greater than P fertiliser application rates and N loading to marine environments was an order of magnitude or two greater than P loading to freshwater. Marine species impacts were higher than freshwater impacts due to strong nutrient delivery; however, strong spatial variations in freshwater fish vulnerability meant freshwater impacts were not necessarily where high loadings were found (e.g. groundnuts versus barley production regions).
At the country level, we identified no significant relationship between crop production quantities and intensities at all metric stages (fertiliser inputs, receptor loadings, aquatic impacts). This suggests there is no benefit from economies of scale to reduce fertiliser impact of eutrophication due to more nutrient management found at all levels of production. Similarly, no strong relationships were identified for the NUE indicator with loading or impact metrics, making the NUE indicator an unsuitable indicator of environmental sustainability from fertiliser use. Our results identify a strong reliance on endpoint indicators that measure nutrient delivery and vulnerability of local receptor environments.
Nutrient transportation was dominated by leaching for N and erosion for P. The underlying concepts to develop the FF suggest that there is an inherent loss of nutrient from fields to receptors. As such, improving nutrient use efficiencies is not sufficient to reduce aquatic impacts. Developing mitigation practice to reduce inherent loss of nutrients through surface runoff, erosion and leaching are important to reduce aquatic impacts caused by fertiliser use. Reducing production or using alternative crops which grow in less vulnerable regions may also reduce the global eutrophication impact. Under the limiting nutrient concept, the FF's spatial persistence showed whilst freshwater persistence of nutrient was greater further inland and upstream of lakes and reservoirs, marine persistence values were high in coastal regions. However, removing the limiting nutrient concept, we identified large loading quantities of N to freshwater, suggesting better nutrient management practices are required on the field, particularly under co-limiting freshwater environments.
Our study provided important improvements to understand the spatially explicit aquatic biodiversity impacts caused by crop commodities beyond the field edge. Our results should be used to enhance and support strategic pollution control policymaking and environmentally sustainable crop-commodity trading and sourcing decisions. The FF model presented here derives spatially explicit FFs for N and P at the global scale from soils and direct emissions using a single methodology at the 5 arcmin resolution. Overall, our results and previous studies suggest assuming diffusive emissions are 10% of direct emissions, is not a sufficiently robust estimate for diffusive emission FFs. Additionally, large FF variations exist between studies, identifying the strong reliance of FF models on robust global nutrient models.
The FFs developed in this study may aid spatially explicit and site-specific LCA nutrient studies from soils. Future research should implement more co-limitation LCIA methodologies, as there are clear limitations to understanding the true biodiversity impact from N and P using the limiting nutrient concept solely. Further research is required to develop FFs for P in marine environments and spatially-explicit XFs and EFs for both N and P within all water environments under co-limitation scenarios. To enhance the FF's practical use, development of reduction factors for nutrient management strategies will support LCA studies and promote good nutrient practices.