Spatially differentiated marine eutrophication method for absolute environmental sustainability assessments

Method assessing absolute marine eutro- phication impacts due to nitrogen emissions Spatially differentiated characterization of marine eutrophication impacts Absolute environmental sustainability assessment of global nutrient emissions Abilitytoidentifycriticalareaswheresafe operating space is exceeded (cid:129) Method shows coherence with observations of hypoxic events. Marineeutrophication and hypoxia causedbyexcess nutrient availability isa growingenvironmentalproblem.Inthis study, we explore marine nitrogen enrichment in the context of Absolute Environmental Sustainability Assessment (AESA), a method combining lifecycle assessment (LCA) with environmental boundariesaiming to compare environ- mental impacts from an activity (product or system) with the safe operating space (SOS) for the activity. Speci ﬁ cally, weaimto increasethespatialresolutionandimprovelifecycle impact assessment(LCIA) modelsformarineeutrophi- cationforuseinAESAs.Byestimatingaproxyofthearealextentofeutrophicationandhypoxiaincoastallargemarine ecosystems (LME), we increased model resolution from 66 LMEs in the original LCIA method to 289 coastal LME subsegments and updated relevant LME parameters tothe newscale(residencetime, bottom water volume, reference O 2 concentration, primary production rates and depths). The new method was tested and validated by comparing the global andspatially differentiated occupationofSOS byglobal nitrogen emissions with observations and itshowed an improved ability to identify critical areas where the SOS is exceeded, in accordance with observations of hypoxic events. Despite limitations such as the estimation of benthic zone volume and low spatial differentiation of environmental boundaries, the method can be used by AESA and LCA practitioners wishing to assess the impact of nitrogen release on marine eutrophication with a higher and more relevant spatial resolution.


H I G H L I G H T S
• Method assessing absolute marine eutrophication impacts due to nitrogen emissions • Spatially differentiated characterization of marine eutrophication impacts • Absolute environmental sustainability assessment of global nutrient emissions • Ability to identify critical areas where safe operating space is exceeded • Method shows coherence with observations of hypoxic events.

G R A P H I C A L A B S T R A C T A B S T R A C T A R T I C L E I N F O Editor: Deyi Hou
Marine eutrophication and hypoxia caused by excess nutrient availability is a growing environmental problem. In this study, we explore marine nitrogen enrichment in the context of Absolute Environmental Sustainability Assessment (AESA), a method combining life cycle assessment (LCA) with environmental boundaries aiming to compare environmental impacts from an activity (product or system) with the safe operating space (SOS) for the activity. Specifically, we aim to increase the spatial resolution and improve life cycle impact assessment (LCIA) models for marine eutrophication for use in AESAs. By estimating a proxy of the areal extent of eutrophication and hypoxia in coastal large marine ecosystems (LME), we increased model resolution from 66 LMEs in the original LCIA method to 289 coastal LME subsegments and updated relevant LME parameters to the new scale (residence time, bottom water volume, reference O 2 concentration, primary production rates and depths). The new method was tested and validated by comparing the global and spatially differentiated occupation of SOS by global nitrogen emissions with observations and it showed an improved ability to identify critical areas where the SOS is exceeded, in accordance with observations of hypoxic events. Despite limitations such as the estimation of benthic zone volume and low spatial differentiation of environmental boundaries, the method can be used by AESA and LCA practitioners wishing to assess the impact of nitrogen release on marine eutrophication with a higher and more relevant spatial resolution.

Introduction
Hypoxia is an escalating environmental problem affecting the world's coastal waters, with severe consequences for marine life (Diaz and Rosenberg, 2008;Vaquer-Sunyer and Duarte, 2008;Zhang et al., 2013;Breitburg et al., 2018;Le Moal et al., 2019;Malone and Newton, 2020). The number of sites where hypoxia has been reported has been increasing exponentially over the last century (Vaquer-Sunyer and Duarte, 2008). This trend is expected to continue, due to the combined effects of eutrophication (i.e. excess nutrient availability, increasing production of organic matter and oxygen demand of coastal systems) and global warming (increasing respiratory oxygen demand and reducing oxygen solubility and ventilation in coastal waters) (Vaquer-Sunyer and Duarte, 2008).
Science of the Total Environment 843 (2022) 156873 In this study, we focus on oxygen depletion caused by excess nitrogen availability, which is primarily a consequence of human activities such as fertilizer application, agricultural nitrogen fixation, sewage discharges and deposition of NO x emissions from fossil-fuel combustion (Sobota et al., 2013;Beusen et al., 2016). We examine marine eutrophication in the context of Absolute Environmental Sustainability Assessment (AESA), which is a method aiming to support decisionmaking in the context of sustainable development. AESA methods compare emissions and associated environmental impacts from an activity (product or system) with an environmental safe operating space (SOS) assigned to the activity. The SOS is delimited by environmental boundaries such as the Planetary Boundaries (Rockström et al., 2009;Steffen et al., 2015). AESA aims to assess whether the environmental impact of an activity remains within the assigned share of SOS (SoSOS). AESAs attempt to cover all potentially relevant impacts on the environment, including marine eutrophication (Bjørn et al., 2020a). An AESA method consists of four main steps: i) identifying environmental boundaries, ii) calculating SOS based on the environmental boundary and reference levels, iii) calculating environmental impacts from emissions, caused by the studied activity and iv) comparing environmental impacts to SOS assigned to the studied activity (Bjørn et al., 2020a).
A substantial limitation of AESA methods is the current application of globally representative average boundaries for non-global processes such as nutrient emissions (Ryberg et al., 2016;Nash et al., 2017) in spite of the known regional variability in ecosystem sensitivities. A recent AESA study developed a spatially resolved method for quantifying impacts of nitrogen emissions on marine eutrophication in relation to the SOS in coastal water (Bjørn et al., 2020b). However, the spatial resolution of the model for coastal water was too coarse in the sense that exceedance of local scale SoSOS (particularly for areas close to shore) might be "hidden" by the inherent averaging of impacts over large marine areas. Moreover, it was assumed that natural oxygen conditions, required to estimate the SOS, correspond to an oxygen saturation of 100 % in the benthic zone. This condition would require full vertical mixing between surface and bottom layer, which is not always the case in coastal areas (e.g., Breitburg et al., 2018).
The AESA method of Bjørn et al. (2020b) builds on the life cycle impact assessment (LCIA) method of Cosme and Hauschild (2017). An LCIA method links emissions from a studied anthropogenic system with impact on the chosen indicator for an environmental impact category (e.g. global warming or marine eutrophication) using characterization factors (CF). A CF represents the environmental fate, exposure and effect of an elementary flow in the environment (Hauschild and Huijbregts, 2015).
The LCIA method of Cosme and Hauschild (2017) is recognized as the best available LCIA method for marine eutrophication as it fills previously existing gaps in the characterization model (Morelli et al., 2018). Moreover, it is recommended by the United Nations' Global Guidance on Environmental Life Cycle Impact Assessment Indicators as the most suitable method for marine eutrophication (UNEP, 2019). However, the method has some limitations (refer to (Morelli et al., 2018, Bjørn et al., 2020b for a full overview). For example, water residence time in the coastal zone relies on literature values and estimations that are uncertain and yet crucial for the estimated size of the FF used in the LCIA method (Cosme and Hauschild, 2017).
In this study, we aim to improve the LCIA models for marine eutrophication for use in AESAs. Taking the AESA method by Bjørn et al. (2020b), building on the LCIA method by Cosme and Hauschild (2017) (hereafter referred to as the "original method"), as a starting point, we 1) modify the model resolution to a scale that better aligns the geographical scale of the environmental problem in question, i.e. eutrophication and hypoxia. 2) update relevant parameters to the new scale (residence time, bottom water volume, reference O 2 concentration, primary production rates and depths), thereby also improving parameters that have been highlighted as uncertain in the original method, and. 3) test and validate the method by comparing the global and spatially differentiated predicted occupation of SOS for global nitrogen emissions with observations of actually occurring oxygen depletion. Apart from its relevance for AESA, a spatially refined method will also advance conventional LCIA modelling in general helping to close a knowledge gap in existing LCIA methods (Morelli et al., 2018).  Table 1 gives an overview of the primary model parameters, the spatial resolution and sources for each of the model components of our method. Fig. 1. The methodological basis for developing and validating the spatially differentiated marine eutrophication AESA method. Note that, in this study, we only look at DIN emitted to coastal water (black arrows in characterization model illustration), while in a full AESA study, one would look at the DIN emitted from human activities ending up in the different environmental compartments (grey dotted arrows in the illustration of the characterization model in the left side of the figure).  Table 1 Overview of the primary parameters of the model components (fate factor (in characterization model), exposure factor (in characterization model) and safe operating space), their spatial resolution and sources, and parameters used for model segmentation. * Aggregated from native resolution to subsegments. 2.1.1. Spatial segmentation Hypoxia resulting directly from human activities is mainly a problem in coastal waters (Bjørn et al., 2020b;US EPA, 2021). Therefore, the first step in creating suitable model resolution was to restrict the 66 large marine ecosystems (LMEs) in the original LCIA model to the 200 m isobath LME i.e. their coastal parts (cLME) (Fig. 2a). Limiting the original LMEs to their coastal parts more than halved the average size of the segments from 1.27 × 10 6 km 2 to 5.6 × 10 5 km 2 . While bringing the model closer to the spatial scale of the problem in question, the size of the cLMEs may still be too large to be relevant for the specific impact category. The relevant size of the model segments was therefore further investigated by examining coastal nitrogen transport patterns.

Method development
The transport of nitrogen and associated concentrations of dissolved oxygen (DO) as well as the spatial extent of eutrophication and hypoxia, depend on various local or regional factors such as advection, air-sea exchange, photosynthesis and chemolithotrophic production (Eldridge and Roelke, 2010). Our model is a regional-scale model with a global coverage building on a relatively simple model structure, compared to more complex regional marine eutrophication models (e.g., Pätsch and Kühn (2008)). Hence, it was not deemed feasible to define the segments considering complex local physical processes and dynamics (see further discussion in Section 3.3). We estimated a relevant size of the model segments by calculating the global average spatial extent for observed large-scale eutrophication events caused by anthropogenic nitrogen emissions. We used high primary production (PP) rate fields (Sathyendranath et al., 2019) associated with large nitrogen emission fluxes (Beusen et al., 2016) to indicate how far an nitrogen emission can propagate in the coastal segment, to determine a proxy for the spatial extent of eutrophication. A high PP rate was defined as an annual average rate of >1563 mg C/m 2 /d (corresponding to a rate two times the global average coastal PP rate), and a large emission as an emission larger than 280 kt/year (corresponding to a point emission contributing with >50 % of the global average emissions to a cLME). Thirteen high PP fields and associated emission points were identified using these criteria, and an arithmetic average area of the high PP fields was calculated. Each cLME segment was then divided by this average area, dictating the number of subsegments necessary for each cLME in order to keep the size of all subsegments similar to size of the average high PP field ( Fig. 2b) (see Supporting Information (SI)-2 for further details).

Calculating marine eutrophication impacts from nitrogen emissions to the coast
As previously mentioned, this work takes the global spatially differentiated marine eutrophication characterization model by  as a basis. A characterization model calculates CFs as the product of a fate factor (FF), exposure factor (XF), and effect factor (EF) respectively, linking the emission of a substance with its impacts in an impact pathway (Fig. 3). The impact pathway describes the cause-effect chain linking environmental drivers (e.g. tomato cultivation) to pressures (e.g. emission of N) and exchanges resulting in changed environmental states (e.g. reduced O 2 concentration in bottom waters) and effects (e.g. loss of marine species richness).
In this study, we consider the cause-effect chain until reduced O 2 concentration in bottom waters and use this as the environmental impact indicator for nitrogen emissions. I.e. we do not include the final effect in terms of possible species loss in the environment, as this would add further uncertainty to the model.
The CFs of the characterization model express the proportional change in environmental impact per change in quantity of environmental interventions (Hauschild and Huijbregts, 2015). In the case of nitrogen emissions and their impact in subsegments of cLMEs, the environmental impact is expressed as shown in Eq. (1):  (Cosme et al., 2015), as there is only PP data for small parts of the cLME available in Sathyendranath et al. (2019) (as most of the cLME is ice). Using the average of these few data points creates heavily overestimated average PP. Fig. 2. Spatial segmentation in the Gulf of Mexico (LME 5) from LME scale to coastal LME scale (cLME) (a) and from cLME to subsegments (b). presence of DIN in the coastal water subsegment to a decrease in oxygen in the benthic zone.
In the elaboration of FF ss , the removal processes by advection (λ a [year −1 ]) and denitrification (λ d [year −1 ]) are considered as first order removal processes. Cosme et al. (2018) suggest that the advection removal coefficient of DIN corresponds to the inverse of the freshwater residence time (Ƭ fw ), i.e. λ a ¼ 1 Ƭfw year ½ and that FF cLME can be written as: Adopting the freshwater residence time concept as an expression for nitrogen removal by advection is inspired by e.g. Dettmann (2001). The residence time concept involves the assumption that all water exchanges (inputs to and outputs from the coastal subsegment) are reflected in this value. See further details on considerations and assumptions on the mass flow balance underlying Equation 3 in SI-1.
In this study, we approximate Ƭ fw to the coastal residence time (CRT) parameter from Liu et al. (2019), which is defined as "the elapsed time since a parcel of source water enters the coastal domain (defined by the 200-m isobath) before it exits to the open ocean while allowing excursions across the coast-ocean boundary" (Liu et al., 2019). The CRT is calculated for each subsegment and is modelled considering a depth and areaweighted average age concentration of an ideal tracer (i.e. a tracer that is not degraded). Hence, only the advective forces influencing the age of the tracer are considered. Liu et al. (2019) provide a high-resolution age concentration dataset (1/8 degree resolution). The 1/8 degree resolved CRT data was aggregated into the 289 cLME subsegments.
The CRT in Liu et al. (2019) reflects the time a water parcel has spent in the entire water column of coastal area. Ideally, it should reflect the residence time in the euphotic zone of the subsegment only, where nitrogen is used in phytoplankton production since this is a perquisite for the oxygen depletion that may cause hypoxia. We assume that the CRT is indeed representative for the residence time in the euphotic zone, considering a water parcel spends more time in the part of the shelf closest to land than the part close to the ocean boundary (Liu et al., 2019). Shallow water depths close to land, in general, result in the entire or most of the water column being within the euphotic zone. In the original method by Cosme and Hauschild (2017), operating with the whole LME segment and not only the coastal part, the above assumption would be too imprecise considering that some LMEs are up to 500 m deep. We believe, however, that this assumption is appropriate at the spatial scale used in our method.
The XF ss kgO2 kgN h i considers plankton growth, the sinking of organic matter and the decrease in dissolved oxygen owing to bacterial degradation of the organic matter as is calculated by Eq. (4): Refer to Table 1 for an overview of the model parameters and (Cosme et al., 2015) for further details. In the original XF model, six of 30 parameters were LME-specific, and the remaining parameters were either expressed at the level of climate zone or global scale. We updated all LME-specific parameters to subsegment scale.

Identification of environmental boundaries and calculation of safe operating space
The environmental boundary forms the basis for calculating the SOS in an AESA assessment. The environmental boundary determines what can be considered a "safe" distance from a threshold or dangerous level of human pressure on the environment (Vea et al., 2020). We applied two environmental boundaries for minimum O 2 concentration in the assessment based on three different boundary approaches (Table 2). Boundary 1 is based on differentiated O 2 levels for five climate zones from Bjørn et al. (2020b) (4.25-7.67 mg/l) (hereafter referred to as boundary approach A), except for oxygen minimum zones (OMZ) where a concentration based on Gallo and Levin (2016) is applied (0.44 mg/l) (hereafter referred to as boundary approach C). Boundary 2 is a global boundary based on Breitburg et al. (2018) (2 mg/l) (hereafter referred to as boundary approach B), except for OMZ where boundary approach C is applied. OMZ are persistent oxygen-depleted areas, and in areas where the OMZ contact the bottom, the benthic species are adapted to lower oxygen concentrations (Diaz et al., 2013). According to Diaz et al. (2013) and Gallo and Levin (2016), cLME #2, 3, 4, 11, 13, 28, 32, 34 and parts of cLME 29 pertain to naturally occurring OMZ (refer to Table S3 (SI-2) for an overview of cLME numbers and names). The boundary approach of these areas (approach C) was defined according to the O 2 tolerance level in OMZ habitats where demersal fish species thrive. We estimated this by assuming that the tolerance level corresponds to the minimum O 2 levels observed in the OMZ where fish species are adapted to lower oxygen levels. This data was retrieved from Gallo and Levin (2016), including only the instances where the minimum depth of the OMZ is <200 m (41 data points in total), in order to exclude open ocean OMZ. This yielded an average O 2 tolerance level of 0.44 mg/l for OMZ subsegments. See SI-4, Table S8 for details.
To ensure transparency and consider the suitability and possible drawbacks of the selected boundary approaches, we analysed them according to the framework in Vea et al. (2020) (Table 2). The framework suggests five components that should be considered when developing or adopting environmental boundaries (refer to Table 2 for an overview and Vea et al. (2020) for a detailed description of the framework and each component).
As seen in Table 2, boundary approach A and B differ in terms of all boundary components, except for the boundary principle, where both boundaries are based on a threshold value. Firstly, they have different overall objectives, where approach A has a more precautionary nature, aiming to prevent harmful effects, while approach B has the more broad objective of avoiding fisheries collapse. The difference is further reflected in their respective response variables and level of accepted impacts and, finally, the value of the control variable (O 2 concentration in bottom waters) which is 2.1-3.8 times higher and more protective for approach A than B. Similar observations are addressed in Vaquer-Sunyer and Duarte (2008), where they showed that empirical hypoxia thresholds vary greatly across marine benthic organisms and that the conventional definition of 2 mg/l is below the empirical sub-lethal and lethal O 2 thresholds for half of the species tested. Breitburg et al. (2018) also argue that, although the 2 mg/l threshold is useful as a generic threshold, a more appropriate approach is to define hypoxia as oxygen levels sufficiently low to affect key or sensitive species. Boundary approach A accounts for spatial differentiation in terms of species sensitivity in different climate zones and is altogether considered more suitable for an AESA study on marine eutrophication impacts than B. However, to ensure feasibility for model validation (see Section 2.2) and to test the influence of the boundary selection on the overall conclusion, boundaries from both studies were considered in this study. Boundary approach C applied for OMZ, differs from A and B considering all boundary components, except the boundary principle. It is, however, most similar to B, where the objective for both approaches B and C is to avoid that fish species disappear. The level of precaution of approach C can be considered as lying between the two other approaches. Boundary approach C is based on the average of 41 fish species, hence more cautious than B, which is based on two species only, but less cautious than A, which is based on the most sensitive species. Despite the different rationale behind boundary C for OMZ, we consider it better to use it in combination with approach A and B than to not account for spatial differences in adaptation of OMZ O 2 levels at all.
For each subsegment, a SOS was calculated as the absolute difference between the boundary (B) (i.e. boundary 1 or 2) and natural reference value (R) (Ryberg et al., 2018;Bjørn et al., 2020b): R represents natural conditions, i.e. the O 2 concentration in the benthic zone before large-scale anthropogenic changes to the nitrogen cycle took place. In cases where B ≥ R, the SOS is zero, as the level of oxygen is already below the boundary. Data on preindustrial O 2 concentration was estimated by subtracting the historical changes in O 2 concentration in bottom water simulated by six CMIP6 models (CanESM5, CNRM-ESM2-1, GFDL-ESM4, MIROC-ES2L, MPI-ESM1-2-HR, UKESM1-0-LL) from present day observations from the WOA, 2018 (Garcia et al., 2019). The historical changes are calculated as the difference between the multi-model average of 2001-2020 and the multi-model average of 1850-1869.
In order to express the SOS in metrics compatible with the characterization model, we converted the SOS from a concentration (in mg/l) to coastal subsegment compartment mass (kg O 2 /subsegment). For the conversion, we estimated total volume of the benthic zone in each subsegment. Where Bjørn et al. (2020b) assumed a generic benthic zone height of Table 2 Overview and analysis according the framework proposed in Vea et al. (2020) of boundaries approaches forming boundary 1 (approach A + C) and boundary 2 (approach B + C). no "safe distance" from threshold and precautionary principle is applied for the final definition of the boundary. However, the method underlying the estimation of the threshold itself (species sensitivity distribution) handles uncertainties and variabilities in data points statistically and tends to include a certain amount of precaution (Belanger et al., 2017). Moreover, the method is based on the most sensitive species and is thus considered precautionary.
The boundary is placed at the threshold value; hence, no "safe distance" from threshold and precautionary principle is applied for the final definition of the boundary. The uncertainty principle underlying the threshold itself is not explicitly stated, however, considering that it is based on the response of only two species and that is it below the empirical sublethal and lethal O 2 thresholds for half of the species tested in Vaquer-Sunyer and Duarte (2008), it is not considered precautionary.
The boundary is placed at the threshold value; hence, no "safe distance" from threshold and precautionary principle is applied for the final definition of the boundary. Considering that it is based on the average O 2 concentration across 41 different fish species, it can be considered more protective than boundary 2 but less than boundary 1. 20 m (as suggested by ), we accounted for spatial differences, relating benthic zone height with total water column depth. For some shallow coastal areas such as the Kattegat region of LME 22 (average depth of 23 m (Oceana, 2014)), assuming a depth of the benthic zone of 20 m is inappropriate. Benthic zone height of 20 m is originally from Doray et al. (2009) and indicates where bentho-pelagic species are most likely to occur in the Bay of Biscay. Ideally, benthic zone height should be defined based on site-specific physical parameters influencing the conditions, i.e., turbidity, temperature and density that make the benthic bottom zone habitable for bentho-pelagic species. For example, Caballero-Alfonso et al. (2015) defined the bottom layer relevant for hypoxia in the Baltic Sea based on pycnocline depth (where the density gradient is the largest in the water column). No study was found with data on physical parameters such as pycnocline depths applicable to estimate site-specific coastal benthic zone heights with global coverage. Therefore, we differentiated the benthic zone height (h) according to total water column depth (d): The factor 0.22 is the average ratio between benthic zone height and total water column estimated from 38 and 33 data points from the Gulf of Mexico (Obenour et al., 2013) and Baltic Sea (Caballero-Alfonso et al., 2015), respectively. Differentiating the original generic benthic zone depth according to total water column depth is the best approach considering the current data availability. However, when such data become available, they can easily be introduced used to further improve the model (see Section 3.3 for further discussion). A list of reference values, benthic zone depths and calculated SOS for all compartments and locations is available in SI-4, Table S4.

Sustainability assessment and method validation
The method was validated by comparing the occupation of SoSOS for environmental impacts on marine eutrophication of global nitrogen emissions with observations of occurrences of hypoxia from Breitburg et al. (2018) and seafloor O 2 levels based on data from WOA 2018 (Garcia et al., 2019). As we in this paper study the nitrogen emissions from all global activities, we assign 100 % of the SOS to these activities. Hence, the SoSOS simply equals the total SOS, and we will refer to it as SOS hereafter. Ideally, the method should predict high occupations of SOS in areas where hypoxia is observed. Environmental impacts (EI) were calculated according to Eq. (1), considering total global nitrogen emission fluxes to the marine compartment (kg/year) from rivers (Beusen et al., 2016(Beusen et al., , data from year 2000 and atmospheric deposition (Ackerman et al., 2019(Ackerman et al., , data from year 2016. The sum of emission fluxes to each cLME subsegment was calculated using the spatial join function in ArcGIS Pro v2.6.2 (ERSI, 2021). Tabulated emission fluxes and associated EI per cLME subsegment can be found in SI-4, Table S4. The occupation of SOS (occ.SOS) was calculated according to Eq. (7) for each cLME subsegment:

Sensitivity and uncertainty assessment
The sensitivity of the updated model parameters (residence time, bottom water volume, reference O 2 concentration and XF) was tested. Note that the individual parameters of the XF were not tested, as this was already tested in Cosme et al. (2015). The sensitivity was calculated as normalized sensitivity coefficients (S coef ) with 25 % perturbation of selected parameters (p) (Ryberg et al., 2015): where p o denotes the original input parameter value, Δp the difference between the default input parameter and the perturbed input parameter, occ. SOS o the original occupation of SOS, and Δocc. SOS o the difference between occ. SOS o and the occupation of SOS calculated for the perturbed input parameter. The sensitivity of a parameter is considered medium if the subsegment average |S coef | ≥ 0.3 and large if the largest |S coef | ≥ 0.5 (Ryberg et al., 2015). We conducted a qualitative uncertainty analysis discussing limitations of the study considering the uncertainty of the model parameters.

Results and discussion
3.1. Spatial differentiation, fate and exposure factors Following the segmentation procedure described in Section 2.1, we divided the global cLMEs into a total of 289 subsegments (SI -2, Fig. S3). This rendered subsegments with an area similar to the global average areal extent of marine eutrophication estimated in this study (in the range of 81,368 km 2 +/− 50 %), except for segment 10, 60, 64 and 65, for which the size of the cLME was already below the range and no subsegmentation was needed. Refer to SI-2 for more details and an overview of number of subsegments per cLME and SI-4, Table S5, for an overview of the area of each subsegment.
Fate and exposure factors were calculated for all 289 cLME subsegments (SI-4, Table S6 and Table S7). The average FF is 0.57 ranging from 0.02 to 2.81 kg N subsegment /kg N emitted × year −1 , similar to the FF in (Cosme et al., 2018) 0.74 (0.03-2.8). The spatial variability between the highest and lowest FF is large relative to the one in Cosme et al. (2018) (181 versus 93) as seen in Table 3.
Low fate factors (below 0.1) are observed in several subsegments located in the West coast of South America (#11.1, 11.2, 13.3 Fig. 4). The low FFs reflect short residence times in these areas. In contrast, we observe high FF (above 1.5) in all subsegments of the Baltic Sea (#23.1, 23.2, 23.3, 23.4, 23.5), North Australian Shelf (#39. 6,39.7,39.8,39.9),Hudson Bay Complex (#63.2,63.6,63.7,63.11,63.12), reflecting higher persistence of nitrogen in these areas (subsegments with the two darkest tones of blue in Fig. 4). Spatial variability between subsegments within cLMEs is on average 2.8, hence relatively small compared to the variability in FF across all subsegments at the global scale (181 as of Table 3). It is largest for cLME 32 with FFs ranging from 0.13 to 1.11 for subsegments #32.8 and 32.2, respectively, with a variability of 14.2. Note that subsegments pertaining to cLMEs #10, 23, 61, 62, 63, 64 are inland seas, island chains, or heavily ice-covered systems for which the estimated residence times, that are the basis for the FF, are uncertain (Liu et al., 2019).
The FFs derived in this study are not directly comparable with those from Cosme et al. (2018), as they represent the persistence, or the length of time a nitrogen load stays in the coastal part of the LME only and not Table 3 Comparison of statistics of FFs and XFs derived in this study and Cosme et al. (2018) and Cosme et al. (2015). FF (kgN subsegment /kg N emitted × year −1 ) XF (kgO 2 /kgN) This study (Cosme et al., 2018) This study (Cosme et al., 2015) Min the whole LME segment as in Cosme et al. (2018). However, if we average the FFs of each subsegment derived in this study over their respective cLMEs, we observe a similar pattern where the FFs for the Baltic Sea, Black Sea, Hudson Bay Complex, Arabian Sea and Yellow Sea (#23, 62, 63, 32, 48) rank high in both studies and Caribbean Sea, Agulhas Current and Southeast U.S. Continental Shelf (#12, 3 and 6) rank low. We also observe discrepancies in particular for Sulu-Celebes Sea, Central Arctic Ocean and Antarctic (#37, 64 and 61), which rank high in Cosme et al. (2018) but low in this study, and Gulf of Thailand, Celtic-Biscay Shelf and Scotian Shelf (#35, 24 and 8) that rank low in Cosme et al. (2018) but high in this study. For #61 and 64 the discrepancy can be due to the high uncertainties entailed to these areas as explained above. For the remaining areas, it can be due to high uncertainties in Cosme et al. (2018), where FFs were estimated based on different studies applying contrasting methods and definitions for estimating residence time.
The average XF in this study is 3.48 ranging from 0.51 to 11.72 (kg O 2 /kg N) with a spatial variability of 23, where Cosme et al. (2015) showed similar XFs (average 5.17 (0.45-15.94)), and a somewhat larger spatial variability (35.4) ( Table 3). As shown in Cosme et al. (2015), PP Pot have the largest sensitivity for the XF model (a sensitivity coefficient of 0.92), and can to a high extent explain the spatial distribution of XFs displayed in Fig , 40.2, 40.4, 39.7, 39.8, 39.9, 45.1, 45.3, 45.4,), East Brazil Shelf (#16.1, 16.2). Comparing with Cosme et al. (2015), XFs for cLME #10, 40, 16, 45, 12 rank among the five lowest in both studies and #23 and 28 rank among the 10 highest. Large discrepancies are observed for e.g. LME #64 (ranking as the lowest in Cosme et al. (2015) and as the 8th highest in this study) and #59 (ranking as number 56 in Cosme et al. (2015) and 9 in this study). For #64 the XF is likely overestimated in our study, as there is only PP data for small parts of the cLME (as most of the cLME is ice). For #59 the reason can be that our study only operates with the coastal part of the LME and PP pot averaged for the whole LME is lower than for the cLME.

Occupation of safe operating space of global emissions
The occupation of SOS for current global nitrogen emissions, applying Boundary 1 and 2, is displayed in Fig. 6. Tabulated results of global nitrogen emissions from rivers and atmospheric deposition to each subsegment can be found in Table S4 (SI-4) along with their occupation of SOS. The SOS is exceeded for 12 subsegments (#5.3, 48.1, 48.4, 14.12, 52.2, 36.1, 62.1, 22.6, 23.2, 47.1, 52.1 and 48.3) and the occupation is high (>60 %) in five subsegments (#52.4, 47.4, 48.5, 27.2 and 17.3), when applying  Boundary 1 (red and orange subsegments in Fig. 6a). When applying Boundary 2, the SOS is exceeded in seven subsegments (#62.3, 47.5, 23.5, 48.1, 48.4, 62.2, 23.4) and the occupation is high (>60 %) in three (#14.12, 36.1 and 23.1) (red and orange subsegments in Fig. 6b). For some subsegments, the natural levels (R) are lower than the boundary (23 and 1 subsegments for Boundary 1 and 2, respectively). For example, subsegment #23.5 (Baltic Sea) has an estimated preindustrial O 2 concentration of 4.24 mg/l, which is lower than Boundary 1 applied for this region (i.e. 7.67 mg/l for subpolar climate zone). For these subsegments, the SOS was set to zero (grey and labeled as "Not applicable" in Fig. 6). Fig. 7 shows a comparison between the occupation of SOS in the original method operating with a spatial scale of 66 LMEs (a) and our method (b), together with maps of seafloor O 2 levels based on data from WOA 2018 (Garcia et al., 2019) (c) and observed hypoxic events (d). Areas where the methods ( Fig. 7a and b) indicate critical low oxygen zones (i.e. high occupation of SOS) in accordance with data from WOA and observed hypoxic events ( Fig. 7c and d) are marked with solid black line in Fig. 7a and b. Both methods capture the observed low O 2 levels and high occupation of SOS in cLME 23, 47 and 48.
However, there are several areas where our method shows exceedance or high occupation of SOS in accordance with WOA and observational data, where the original method does not. Areas where the original model fails to capture areas of low oxygen indicated in Fig. 7c and d are marked with a dotted line in Fig. 7a. For instance, around the Arabian sea (#32.8), Bay of Bengal (#34.2), West Coast of America (#4 and 11.2 marked as "not applicable" in our method, as the preindustrial levels were already below the boundary), the Red Sea (#62), Northeast coast of Africa (#31 and 32.5), and Northern part of the Gulf for Mexico (#5.3). This indicates that our method is better at predicting low O 2 levels because its resolution is closer to the spatial scale of eutrophication based oxygen declines. For example, about 70 % of the river flux to the Gulf of Mexico (cLME #5) comes from the Mississippi River located in the northern part of the Gulf (subsegment #5.3). According to observations in Fig. 7d, the majority of hypoxic events are located here. The potential contribution to hypoxia from this point emission is not captured in the original method, as the impact from the emission is diluted over the volume of the whole LME. This indicates that the environmental mechanisms and impacts caused by a point source occur in areas smaller than the full LME, which was the basis for the original method. In Gulf of California and the coastline of Peru (subsegment #4 and #13.3), our model also fails to capture the areas of low oxygen indicated in Fig. 7c and d (stippled line in Fig. 7a and  b). However, these areas are OMZ where we used a boundary lower than 2 mg O 2 /l (0.4 mg O 2 /l). Hence, the maps of Fig. 7c and d indicating critical areas below 2 mg O 2 /l are not comparable with the models for these areas. Refer to Section 2.1 for explanation and rationale for a differentiated boundary for OMZ.
Finally, there are some areas where there is inconsistency between Fig. 7c and d, where one map indicates critical O 2 zones and the other does not. This may be due to uncertainties underlying these maps such as lack of observational data in developing countries (Breitburg et al., 2018), which may be the case for the Southwest coast of Africa (subsegment #29.2). Fig. 7d includes observed seasonal hypoxic events where Fig. 7c shows the yearly average O 2 concentration only, which can explain the discrepancies at e.g. the Southwest and Southeastern coast of Australia (#43.1, 42.2 and 42.3). Moreover, some areas may be subject low nitrogen emissions not leading to observed hypoxic events or low occupation of SOS, despite that these areas have naturally low oxygen concentrations. For example, Oyashio Current and Aleutian Islands (subsegments #51 and 65) show low occupation of the SOS (below 1 %) in Fig. 7a and b, no observed hypoxic events in Fig. 7d, but low oxygen concentrations in Fig. 7c. In these areas, our data show low natural oxygen levels (3.8 mg/l for #51 and 2.8 mg/l for #65), hence in accordance with Fig. 7c. However the nitrogen loads to these cLMEs are also very small (only 5.5 % and 0.2 %) of the average cLME nitrogen load), hence they are not shown as critical oxygen areas with high occupation of SOS in our model. Areas where the models do not indicate high occupation of SOS and low oxygen areas, shown in either Fig. 7c or 7d, are marked with a dotted line in Fig. 7a and b.

Method limitations and future research needs
We found that the results have a high sensitivity for all tested parameters (Table 4), with highest sensitivity to the O 2 reference level followed by XF, volume and CRT.
The O 2 reference was estimated based on observations of current O 2 levels from WOA 2018 (Garcia et al., 2019) and historical changes in O 2 concentration in bottom oceans simulated by six CMIP6 models. The WOA data is relatively coarse (1 degree) and is not resolved on fine-scale coastal features. Moreover, there is uncertainty regarding historical O 2 level changes, simulated by the six numerical models. Numerical models tend to simulate a decline in the total global ocean oxygen equal to only about half of the most recent observation-based estimates. Moreover, they predict different spatial patterns of oxygen changes, in particular, in the tropical thermocline (Breitburg et al., 2018). However, as we apply an average of simulations of six models, the uncertainty is considered acceptable.
The most sensitive parameter of the XF model was in Cosme et al. (2015) shown to be the PP pot . The PP rates we used were in the original source estimated by translating remotely sensed surface ocean optical properties to PP. Translating these observations into PP is not straight forward (Gregg and Rousseaux, 2019;Sathyendranath et al., 2019). For example, an important factor for estimating phytoplankton is photo-adaptation, which is difficult to estimate from ocean color satellite observations (Gregg and Rousseaux, 2019). Moreover, the resolution is relatively low to represent detailed dynamics in coastal areas (1 degree). Future method improvements should aim to obtain data on PP with less uncertainty and higher resolution. For example, Mattei and Scardi (2021) updated and improved a global dataset on marine phytoplankton primary production, handling missing data and adding variables related to primary production, which were not present in the original datasets. When available, the method presented in this study could be updated using the new data.
The volume applied to estimate the SOS was based on the area of the subsegments and benthic zone depth. Differentiating the benthic zone depth according to total water column depths is the best approach considering the current data availability. However, this is a simplified and generic approach not considering spatial differences in the ratio between the benthic zone and total water column depth. When data on physical parameters such as the pycnocline depths applicable to estimate site-specific coastal benthic zone heights with global coverage is available, future improvements of the method should focus on this. Despite having the lowest average sensitivity of the variables tested, the sensitivity of the CRT is also high. Even though the parameter has uncertainties entailed as discussed in Section 3.1, it is considered the currently best available data with a global coverage and high resolution.
Although not explicitly tested in the sensitivity analysis, we showed in Section 3.2, that the selection of boundaries delimiting the SOS has a high influence of the results and, for some subsegments, natural levels are lower than the boundary delimiting the SOS, making the method not applicable for these areas. We accounted for higher tolerance of low O 2 levels in  OMZ, however, there are also other areas where O 2 levels are naturally low due to, for example, high residence times. The spatial differentiation of the boundaries should increase in order to account for different O 2 tolerance in such areas as well. Moreover, boundary approach C used for OMZ should be revised and be based on the same principles as the boundary approaches it was combined with. For example, for boundary approach A from Bjørn et al. (2020b) it should be based on LOEC concentrations. Acknowledging that there are spatial differences in terms of the areal extent of eutrophication, dictated by variations in the physical transport mechanisms, we used a global average areal extent across large-scale events for the segmentation. Obtaining spatially differentiated segments describing areal extent of nitrogen transport from a point source would require high resolution spatial data describing local effects such as transport flows and directions influencing freshwater plume structures. These local effects associated with regions of freshwater influence from rivers and local point sources could be addressed with high resolution CRT fields. However, CRTs with sufficiently high resolution to describe these local effects (higher than applied in this study, i.e., 1/8 degree) is currently not available with global coverage for coastal segments.
We have shown that the resolution of our method allows for identifying critical areas where the SOS is exceeded, in accordance with observations and hence caters to the modelling of environmental mechanisms involved in the marine eutrophication impact pathway. However, for determining the optimal resolution, one must consider the purpose of the study. For an LCA or AESA study, we argue that our model is sufficiently resolved since a higher resolution would be less compatible with the resolution of typical life cycle inventory (LCI) data used in such studies. The geographical resolution of LCI data e.g. from the ecoinvent database (Wernet et al., 2016), is available at state, country or continent or global level. Further research should test our method and its resolution in a case study considering whether and how CFs could be aggregated to match LCI data with a coarser spatial resolution.
Our study shows that substantial improvements in the impact assessment of marine eutrophication can be obtained by increasing the spatial resolution, and we hope that this can inspire to improve LCIA methods of other impact categories as well by updating their spatial resolution to a more relevant one. We also encourage the LCIA and AESA community to use our method validation as inspiration to a more general validation of LCIA methods with measurements, which is often lacking, to improve the overall credibility of the impact assessment in LCA.

Conclusions
This study substantially improves LCIA and AESA of marine eutrophication by modelling impacts at relevant spatial scales and by refining the LCIA model with updated and spatially representative data. We increased the model resolution from 66 LMEs to 289 coastal LME subsegments representing the geographical scale of eutrophication and hypoxia, and updated relevant parameters to the new scale (residence time, bottom water volume, reference O 2 concentration, primary production rates and depths).
A model validation comparing model outputs with observations of hypoxic events showed that our method is able to better predict areas with critical oxygen concentrations. Despite its limitations, the method can be used by AESA and LCA practitioners that wish to assess impact on marine eutrophication with a higher and more relevant spatial resolution. We believe that the improved estimates of marine eutrophication impacts from human induced nitrogen emissions can support more informed and better decisions for managing nitrogen and reducing emissions to within a sustainable level.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.