Riverine nitrogen supply to the global ocean and its limited impact on global marine primary production: a feedback study using an Earth system model

A common notion is that negative feedbacks stabilize the natural marine nitrogen inventory. Recent modeling studies have shown, however, some potential for localized positive feedbacks leading to substantial nitrogen losses in regions where nitrogen fixation and denitrification occur in proximity to each other. Here we include dissolved nitrogen from river discharge in a global 3-D ocean biogeochemistry model and study the effects on near-coastal and remote-open-ocean biogeochemistry. We find that at a steady state the biogeochemical feedbacks in the marine nitrogen cycle, nitrogen input from biological N₂ fixation, and nitrogen loss via denitrification mostly compensate for the imposed yearly addition of 22.8 to 45.6 Tg of riverine nitrogen and limit the impact on global marine productivity to < 2 %. Global experiments that regionally isolate river nutrient input show that the sign and strength of the feedbacks depend on the location of the river discharge and the oxygen status of the receiving marine environment. Marine productivity generally increases in proximity to the nitrogen input, but we also find a decline in productivity in the modeled Bay of Bengal and near the mouth of the Amazon River. While most of the changes are located in shelf and near-coastal oceans, nitrogen supply from the rivers can impact the open ocean, due to feedbacks or knock-on effects.

1 Introduction

Nitrogen plays a key role in marine biogeochemistry in coastal and open oceans as it is one of the major limiting nutrients for algal photosynthesis. Variations in the oceanic fixed-nitrogen (N) inventory are known to have driven marine productivity changes contributing to atmospheric CO₂ variations in Earth's history (Falkowski, 1997).

Although several studies have questioned the stability of the global N budget (Codispoti et al., 2001; Gruber and Sarmiento, 1997; Codispoti, 1995), the present marine N inventory is generally considered to be in a steady state (Deutsch et al., 2007; Altabet, 2006; Gruber, 2004; Tyrrell, 1999; Redfield et al., 1963). Oceanic fixed-nitrogen concentrations are mainly controlled by the balance between denitrification and N fixation, creating negative feedbacks which damp the often strong perturbations of the marine N content with respect to the more slowly overturning phosphorus (P) inventory (Somes et al., 2013; Deutsch et al., 2007; Gruber, 2004; Ruttenberg, 2003).

In regions where fixed N is sparse, organisms that fix atmospheric nitrogen (N₂) dissolved in seawater, commonly called diazotrophs, can compensate for N deficits (Deutsch et al., 2007, 2001; Capone et al., 1997). However, as diazotrophs consume phosphate while adding fixed N to their environment, they are generally considered to regulate their own population. Previous studies have shown that at least the most commonly cultured diazotrophs, especially Trichodesmium, have low growth rates relative to many non-fixing phytoplankton (e.g., Capone et al., 1997). Indeed, while these organisms are able to fix N₂ when reactive nitrogen is scarce, this turns into an disadvantage in regions where N is more abundant, because N₂ fixation requires more energy (Tyrrell, 1999). According to common conceptual models of controls on nitrogen fixation, the slowly growing diazotrophs are then outcompeted by non-fixing phytoplankton if enough P and other nutrients are present (Tyrrell, 1999). Note that part of our knowledge of N₂ fixation and most modern model concepts are still based on the original limited assumptions based on a few species, especially Trichodesmium. Also, emphasis has traditionally been put on bottom-up controls, despite accumulating evidence that top-down controls such as selective grazing on dominant species may have substantial impacts on the distribution of diazotrophs and nitrogen fixation (Landolfi et al., 2021).

Denitrification is a metabolic process in which nitrate (NO₃) replaces oxygen (O₂) as the terminal electron acceptor for respiration and is reduced to N₂, which is not bioavailable for most marine organisms, except diazotrophs (Gruber, 2004; Deutsch et al., 2001). Denitrification represents the main sink for fixed N in the ocean and occurs both in marine sediments and in the water columns under suboxic conditions (Gruber, 2004; Codispoti et al., 2001). But denitrification limits itself by reducing the concentrations of fixed N at the surface, which in turn limits the growth of phytoplankton and the heterotrophic O₂ consumption during organic matter remineralization, eventually making NO₃ less competitive as an electron acceptor where O₂ remains available (Landolfi et al., 2013; Gruber, 2004). These two processes, N₂ fixation and denitrification, both contribute to regulating the global marine N budget.

Besides the fixation of atmospheric N₂, rivers are also a major source of N to the coastal and the open ocean. Rivers are estimated to add 36–60 Tg N yr⁻¹ to the coastal waters (Beusen et al., 2016; Mayorga et al., 2010; Seitzinger et al., 2005). These N inputs are regionally highly diverse and range over several orders of magnitude (Meybeck et al., 2006).

Although riverine N is not the main source of N to the marine environments, it can become a key player, as it is directly influenced by human activities. Seitzinger et al. (2010), for example, estimated that global nitrogen export by rivers to the coastal waters increased by 17.7 % from 1970 to 2000. Nitrogen is known to impact coastal marine biology and biogeochemistry, leading for example to eutrophication, algal blooms, or hypoxia (e.g., Seitzinger et al., 2010; Billen and Garnier, 2007; Smith et al., 2003). Previous studies have shown that nutrient input from land also has consequences for seawater composition and through this impacts biogeochemical processes in the open ocean farther away from the coasts (e.g., Barron and Duarte, 2015; Bauer et al., 2013; Bernard et al., 2011; Jahnke, 2010).

At this place it is appropriate to include some remarks about the anthropogenic perturbation of the N cycle. In the Anthropocene, human activities have led to increased inputs of fixed nitrogen from the atmosphere and through different sources of runoff from land (e.g., Somes et al., 2016; Kim et al., 2014; Lamarque et al., 2013). At the same time, warming and deoxygenation can lead to increased N loss (Oschlies et al., 2019). The question has therefore been raised of if the global N budget could still be considered to be in a steady state. While these combined effects on the N budget are still very uncertain, some studies suggest that imbalances could be limited due to internal feedbacks of the N cycle (Landolfi et al., 2017; Somes et al., 2016; Krishnamurthy et al., 2006).

As global measurements of N concentrations and fluxes are difficult, models are often used to study the marine N cycle and its feedbacks. However, global biogeochemical ocean models often still omit riverine nutrient input to the ocean or represent it in a very simplified form (Séférian et al., 2020). Giraud et al. (2008), for example, tested the sensitivity of the global ocean biogeochemistry to coastal nutrient fluxes in a global ocean biogeochemistry model by introducing nutrients in different scenarios into the coastal grid boxes. They found that excess nutrients in the coastal ocean could impact the biological activity not only locally but also in the open ocean and that the effect depended more on the ratio between these nutrients and iron and silicate than on the actual quantities. Nevertheless, the study by Giraud et al. (2008) was an idealized experiment without observed nutrient fluxes and with a relatively simple representation of the ecosystem dynamics, where total nitrogen nutrient and phosphate were linked by the Redfield ratio and indifferently represented by one model variable. Da Cunha et al. (2007) used an ocean biogeochemistry model to analyze the impact of river nutrient fluxes (N, Si, Fe, and carbon) on the global and coastal ocean primary production but concentrated on a short time period of a few decades, likely not long enough to study the feedbacks of the N cycle in the open ocean considering that the mean residence time of fixed nitrogen in the ocean has been estimated to be a few thousand years (Gruber, 2004).

More recently, Lacroix et al. (2020) implemented estimated riverine nutrient loads in a global ocean model to analyze their implications for global oceanic nutrient concentrations, primary production, and CO₂ fluxes. Their focus was on pre-industrial nutrient input from rivers, estimated as a function of precipitation, surface runoff, and temperature. N was calculated from the simulated P using a fixed N : P ratio, but Lacroix et al. (2020) did not analyze the N-cycle feedbacks.

For our study, we used the Earth system climate model of intermediate complexity of the University of Victoria (UVic), version 2.9 (Eby et al., 2009; Weaver et al., 2001). Earth system models of intermediate complexity (EMICs) have been developed to fill the gap between more abstract conceptual models and comprehensive global and Earth system models (ESMs) (Claussen et al., 2002). EMICs allow the integration of a large number of processes, more than conceptual models do, often using a coarser resolution and simplifying assumptions, e.g., describing the atmospheric circulation compared to ESMs, and thus substantially reducing computational costs. UVic has been developed as a tool that helps to understand processes and feedbacks operating within the climate system on decadal and longer timescales (Weaver et al., 2001).

Atmospheric deposition is known to be another important source of N to the ocean. Although it is estimated to add nitrogen at the same magnitude as the rivers and will also become more important with increasing anthropogenic activities (e.g., Tyrell, 1999; Cornell et al., 1995), it will not be considered in this study. Previously, Landolfi et al. (2017) and Somes et al. (2016) used UVic to study the response of the marine N cycle to idealized atmospheric N deposition and its impact on marine productivity. While Somes et al. (2016) performed a series of idealized sensitivity experiments to evaluate the spatial and temporal scales of N-cycle feedbacks, Landolfi et al. (2017) used an atmospheric N deposition forcing reconstructed using the multimodel mean of the Atmospheric Chemistry and Climate Model Intercomparison Project (Lamarque et al., 2013). Both found that N-cycle feedbacks stabilize the model's marine N inventory and limit global changes to the marine N cycle and productivity. But none of these studies included riverine N supply.

In order to disentangle the effects of the different sources of N, we are using UVic without atmospheric N deposition to focus on the marine biogeochemical response to riverine N inputs to the coastal ocean. To do this we make use of modeled estimates of riverine dissolved inorganic nitrogen (DIN) export from watersheds (Mayorga et al., 2010). While atmospheric deposition is more spread out over the whole ocean, river export of dissolved N reaches the ocean as point sources at different locations and in different concentrations. The global amount of N added to the ocean is comparable between our river-supply study and the atmospheric-supply ones by Somes et al. (2016) and Landolfi et al. (2017). Nevertheless, we hypothesize that the response of the marine ecosystem differs with highly concentrated nutrient injections associated with individual rivers.

In order to test the N-cycle mechanisms and feedbacks found and described before, we set up an experiment where we simulate differential riverine nitrogen supply to the coastal oceans. Although riverine nitrogen supply is highly influenced by anthropogenic activities, our focus is on the natural nitrogen cycle. Analogously to earlier studies (Lacroix et al., 2020; Da Cunha et al., 2007), we first evaluate the global N inventory and marine primary production after sustained addition of riverine DIN. In a second step, we additionally performed a series of experiments where we studied the responses of the ocean to riverine nutrient supply to individual regions in order to find out if N impacts the global ocean differently, depending on the region where river supply takes place.

2 Model description and experimental design

Nutrients from the Global Nutrient Export from WaterSheds 2 (NEWS 2) model (Mayorga et al., 2010) are added to the University of Victoria Earth System Climate Model (UVic) 2.9 (Keller et al., 2012; Eby et al., 2009; Weaver et al., 2001). The model is outlined, before describing the NEWS 2 data set and our experimental design below.

2.1 The Earth system model UVic 2.9

UVic (Weaver et al., 2001) version 2.9 (Keller et al., 2012; Eby et al., 2009) is an Earth system model of intermediate complexity (Claussen et al., 2002). It consists of a three-dimensional ($\mathrm{1.8}{}^{\circ }\times \mathrm{3.6}{}^{\circ }$, 19 levels) general circulation model of the ocean: a two-dimensional, single-layer energy–moisture balance atmospheric model; a dynamic–thermodynamic sea ice model; and a terrestrial vegetation model.

The atmospheric component dynamically calculates heat and water fluxes between the atmosphere and the ocean and land and sea ice and is forced by monthly climatological winds prescribed by National Centers for Environmental Predictions/National Center for Atmospheric Research (NCEP/NCAR). The 19 vertical levels of the oceanic component, Modular Ocean Model 2 (MOM2), are 50 m thick near the surface and up to 500 m in the deep ocean. The oceanic physical settings are the same as in Keller et al. (2012). The marine ecosystem module of UVic is based on Keller et al. (2012) with updates of some of the equation parameters as noted in Partanen et al. (2016), where a small error in the code was corrected. Seven prognostic variables are embedded within the ocean circulation: two phytoplankton classes (nitrogen-fixing diazotrophs P_D and other phytoplankton P_O), zooplankton (Z), sinking particulate detritus (D), nitrate (NO₃), phosphate (PO₄), and oxygen (O₂) (Fig. 1). NO₃ and PO₄ are linked through exchanges with the biological variables by constant (Redfield) stoichiometry of organic matter. All biological variables as well as the detritus are expressed in terms of nitrogen (mmol N m⁻³), using Redfield stoichiometry to calculate carbon and P. Since diazotrophs can fix nitrogen gas dissolved in seawater, they are not limited by NO₃ or by a maximum NO₃ concentration, while the growth of other phytoplankton is limited by NO₃ and PO₄ (note that both are additionally limited by iron, light, and temperature). The explicit integration of diazotrophs permitting the computation of nitrogen fixation is not given in all ocean models but makes UVic a good choice to study nitrogen cycle feedbacks. The maximum potential growth rate of diazotrophs is based not only on temperature as in most models but also on dissolved iron, which is necessary, for example, for photosynthesis or the reduction of nitrate to ammonium (Keller et al., 2012; Galbraith et al., 2010). Keller et al. (2012) found that the observational estimates were within the range of global nitrate fixation rates from estimations and that the patterns of N₂ fixation from the new model were mostly consistent with the relatively sparse available observations (Sohm et al., 2011). See Keller et al. (2012) for a full description and evaluation of simulated marine biogeochemistry.

Ecosystem model

Figure 1 Ecosystem model schematics for the NPZD model with the prognostic variables (in square boxes) and the fluxes of material between them, indicated by arrows. The prognostic variables include two nutrients, nitrate (NO₃) and phosphate (PO₄), two phytoplankton (nitrogen fixers P_D and other phytoplankton P_O) as well as zooplankton (Z), sinking detritus (D), and dissolved oxygen (O₂). Nitrate (NO₃) and phosphate (PO₄) are linked through exchanges with the biological variables by constant (Redfield) stoichiometry of organic matter. See further details in the text. Figure updated from Keller et al. (2012).

In the global ocean, fixed N is regulated by the major input fluxes, N₂ fixation and riverine input, and the major removal flux, denitrification (here implicitly including anammox). Benthic denitrification in particular is believed to be the major sink for fixed N (Voss et al., 2013; Galloway et al., 2004). It is included here through empirical transfer functions derived from benthic flux measurements (Bohlen et al., 2012). The functions are based on dynamic vertically integrated sediment models and estimate denitrification from the rain rate of particulate organic carbon to the seafloor and bottom-water O₂ and NO₃ concentrations. Like Somes and Oschlies (2015) and Somes et al. (2013) we use a subgrid bathymetry scheme for shallow continental shelves and other topographical features that are too fine to be resolved on the coarse UVic grid in order to better resolve particulate organic matter sinking and remineralization at the seafloor. For each cell near the coast this scheme calculates the seafloor area within the cell at a higher resolution following Somes et al. (2010b).

2.2 Including riverine nutrient supply to the UVic ocean

2.2.1 Global Nutrient Export from WaterSheds 2 – NEWS 2

Riverine N added to UVic has been generated by a global, spatially explicit model of nutrient exports by rivers. NEWS 2 (Mayorga et al., 2010) is the second version of a system of sub-models, which estimate the present-day annual export yield for each river basin (kg N km⁻² yr⁻¹) at the river mouths for dissolved and particulate forms of organic and inorganic N and P, as well as for dissolved organic and particulate carbon. See Mayorga et al. (2010) for more details on the model configuration. Each sub-model predicts the river export of a nutrient element for the base year 2000. This export is calculated as a function of natural and anthropogenic biogeophysical properties of each of the 5761 exoreic basins considered (Seitzinger et al., 2005). NEWS DIN includes DIN from sewage point sources, as well as N from diffuse sources, mobilized from watershed soils and sediments (Dumont et al., 2005). Despite uncertainties and errors, NEWS DIN predicts 54 %–78 % of the variability in the DIN export yield (kg N km⁻² yr⁻¹) and 72 %–83 % of the DIN export load (kg N per basin per year) of the validation data set used by Dumont et al. (2005). Note that NEWS 2 excludes runoff from the Antarctic continent.

2.2.2 NEWS DIN for UVic

To estimate the total export per river mouth, we multiplied the yields (kg N km⁻² yr⁻¹) of DIN and dissolved organic nitrogen (DON) by the respective basin area (in km²). Data from the NEWS 2 models have been interpolated on the coarser UVic grid, and the total exports per river basins have been added to the nearest discharge points, as not every river mouth from NEWS 2 has its equivalent discharge point in UVic (Fig. 2). Because there can be strong seasonal variations in nutrient fluxes and fluvial nutrient imports can have different effects on the biogeochemistry of a coastal ecosystem depending on the timing of the fluxes (Eisele and Kerimoglu, 2015; Holmes et al., 2012; Townsend-Small et al., 2011), we used the seasonally cycling climatology of freshwater runoff from UVic to estimate seasonal variations in N supply. Although freshwater discharge and riverine nutrient export are not always correlated, the discharge has an important impact on the nutrient loads of rivers (e.g., Lu et al., 2011, 2009; Sigleo and Frick, 2007; DeMaster and Pope, 1996). Here, we assumed a constant seasonal cycle in runoff and that nitrogen concentrations in the discharged river water are constant throughout the seasonal cycle. We then distributed the annual load over the months, weighted by the fraction of monthly freshwater discharge.

Riverine phosphorus is not added in this experiment, and we assume a fixed marine P inventory, like in most previous studies with UVic. The inclusion of a dynamic P cycle (like in Niemeyer et al., 2017; Kemena et al., 2019) with riverine P supply from NEWS 2 will be the subject of a follow-up study.

Figure 2 DIN export yield for each discharge point in mmol m⁻² yr⁻¹ from the NEWS 2 data set interpolated on the UVic grid for January (a), April (b), July (c), and October (d).

2.3 Experimental design

To analyze the effect of riverine nutrient export in the UVic model, three experiments have been performed: NEWS, DIN+DON, and 2×DIN (Table 1). All simulations were run for 10 000 years with benthic denitrification and subgrid bathymetry, starting from an already-spun-up steady state with the standard model version (i.e., with no riverine nutrient input) with pre-industrial conditions for insolation and a fixed atmospheric CO₂ concentration of 283 ppm (Keller et al., 2012). In NEWS only DIN from NEWS 2 was added to the discharge points in UVic. Here, we evaluate how the model and especially the nitrogen cycle react to the riverine nutrient input. In DIN+DON we added DIN and DON from NEWS 2. In 2×DIN twice the yield of DIN from NEWS 2 has been added. These simulations with increased N supply are scaling experiments to test the mechanisms and feedbacks described for example by Landolfi et al. (2017) and Somes et al. (2016). For comparison, a control simulation has been run for 10 000 years without riverine DIN supply (referred to as CTR). Globally, NEWS 2 predicts a riverine N supply of 22.8 Tg N yr⁻¹ for DIN and 11.8 Tg N yr⁻¹ for DON. Both enter the biogeochemical model as NO₃ fluxes in mol N m⁻² s⁻¹, thereby implicitly assuming that all DON is bioavailable or rapidly turned over to DIN. The marine ecosystem dynamics as well as the biogeochemical cycles of the model run have been evaluated in previous studies under the standard boundary conditions, without riverine nutrients (e.g., Somes et al., 2013; Keller et al., 2012; Somes et al., 2010b; Schmittner et al., 2008, 2005). We therefore concentrate on the evaluation of the response of the marine biogeochemical model to the new model component of riverine nutrient discharge. The global ocean biology reacts to the new N-cycle components in the first 3000 to 4000 years of the simulations. After this first phase, the N budget slowly equilibrates towards a steady state (see figure in the Supplement). For the evaluation of the resulting ocean biogeochemistry, we analyze in the following the mean of the last 100 years of each 10 000-year-long simulation.

Table 1 UVic simulations and global annual nutrient flux from river discharge from Global Nutrient Export from WaterSheds 2 (Mayorga et al., 2010).

3 Results and discussions

3.1 Nitrogen

3.1.1 Global nitrate distribution

In comparison with observational data of the World Ocean Atlas (Garcia et al., 2019), the model simulates the general structure of the profiles fairly well but underestimates the observed NO₃ in the water column in each ocean basin by 3 to 4 mmol m⁻³, especially omitting the midwater maximum at around 1000 m (Fig. 3). Global average NO₃ concentrations only vary a little between the simulations (from 22.19 mmol m⁻³ in CTR to 22.48 mmol m⁻³ in NEWS and 22.84 mmol m⁻³ in 2×DIN), and the differences from CTR correspond to +1.1 %, +1.8 %, and +2.5 % of the total observed inventory for NEWS, DIN+DON, and 2×DIN, respectively. Nevertheless, in all three simulations (NEWS, DIN+DON, and 2×DIN), NO₃ concentrations are globally higher compared to CTR (Figs. 5, 6). The absolute error between model and observations decreases with higher riverine N supply (top right panel in Fig. 3). At the surface, the global ocean NO₃ distribution patterns are very similar between the model and the observations, as well as between the control (CTR) and the NEWS simulation (Fig. 4).

Figure 3 Profiles of global averaged NO₃ in mmol m⁻³ from UVic simulations with and without riverine DIN export from NEWS data sets and observations with the WOA. (a) Global ocean average of NO₃. (b) Global profiles of misfit in NO₃ compared to the observations. (c–e) NO₃ profiles for the three main ocean basins: the Atlantic Ocean (ATL), Pacific Ocean (PAC), and Indian Ocean (IND).

Figure 4 Surface NO₃ concentrations for simulations UVic NEWS (a) and UVic control (b) and from observations from the World Ocean Atlas (Garcia et al., 2019) (c).

Figure 5 Difference in zonal mean ocean concentrations of NO₃ between the UVic simulations with riverine DIN export and the control simulation. (a, b, c) Difference between NEWS and CTR; (d, e, f) difference between DIN+DON and CTR; (g, h, i) difference between 2×DIN and CTR. The columns show the zonal mean of the Atlantic Ocean (a, d, g), the Pacific Ocean (b, e, h), and the Indian Ocean (c, f, i) basins. The difference in zonal averaged NO₃ concentrations is higher than the color bar maximum at the surface in the northern Indian Ocean basin with a maximum for 2×DIN at 7.2 mmol N m⁻³. Note that the three ocean basins have different sizes in terms of latitudes but for layout reasons the panels have the same dimensions.

The supply of riverine NO₃ affects the ocean nutrient concentration not only locally near the river mouths but also in regions far away from the coasts. Surface NO₃ concentrations increase with higher river supply in the coastal regions and in the higher latitudes. The globally highest increase in NO₃ can be found in the 2×DIN experiment (see also Table 2), and the NO₃ increase is higher in the deeper ocean than at the surface. Interestingly, the increase in the oceanic N inventory is more than twice as high in UVic 2×DIN compared to UVic NEWS, indicating non-linear feedbacks.

While higher NO₃ concentrations due to riverine input are not entirely surprising, some regions present lower concentrations compared to CTR. In all simulations NO₃ is slightly lower at the surface in large parts of the tropical and subtropical oceans. At 850 m depth, the ocean loses NO₃ upon the addition of riverine N in our simulations in low-oxygen regions where denitrification occurs, such as the Benguela Current system, the Bay of Bengal, and the eastern equatorial Pacific near the coast of Central America (Fig. 6).

Figure 6 Difference in global distribution of NO₃ concentrations between experimental simulation and the control simulation in mmol N m⁻³. The different simulations are NEWS, DIN+DON, and 2×DIN. The differences are shown at three different depths: surface (top), 850 m (middle), and 1100 m (bottom).

Table 2 Amount of additional nitrogen in teragrams of N at the end of each simulation compared to CTR.

In major parts of the Atlantic Ocean and Pacific Ocean basins, NO₃ concentrations are higher in NEWS than in CTR (Fig. 5). NO₃ is particularly elevated in the upper 2000 m in the North Atlantic Ocean (up to 2 mmol N m⁻³ in 2×DIN, corresponding to +18 %) and upper 1000 m in the North Pacific, but the difference between the simulations is positive across the whole of both basins, indicating that a substantial part of the additional riverine N is exported into the open and deeper ocean. At the surface of the tropical and subtropical oceans, however, NO₃ concentrations are lower by a maximum of 0.9 mmol N m⁻³ in the UVic NEWS experiments compared to CTR (Fig.  6).

The Indian Ocean basin comprises the Arabian Sea and the Bay of Bengal. Zonally averaged NO₃ concentrations reflect essentially the behavior of the Bay of Bengal, where the rivers of the Ganges Delta supply high amounts of nutrients to the northern basin (Figs. 5, 7). Here, the model simulates high NO₃ concentrations at the surface (in the north, several hundred percent higher in NEWS when compared to CTR). But in the deeper northern Indian Ocean basin down to approximately 2000 m, NO₃ concentrations are significantly lower in NEWS than in CTR. Considering the zonal average of the Indian Ocean, NO₃ concentrations are lower by 0.7 to 0.9 mmol N m⁻³ and even more if only the zonal average of the Bay of Bengal is considered.

Part of the global NO₃ patterns can be explained by the interaction of ocean circulation and biology. N is transported into the interior ocean via circulation and also accumulates due to the biological pump. But these processes do not explain the loss in N in the subtropical surface oceans or the Bay of Bengal. The same applies to the total amounts of N. Despite the continuous supply from the rivers, the additional NO₃ in NEWS, DIN+DON, and 2×DIN compared to CTR amounts to an increase of only 1.1 %, 1.8 %, and 2.5 %, respectively (Table 2). What limits the increase in the global oceanic N inventory is the combination of the N-cycle processes denitrification and N₂ fixation.

Figure 7 Difference in zonal average of NO₃ and O₂ in the basins of the Arabian Sea (a, b, e, f) and the Bay of Bengal (c, d, g, h). The panels in the left column show the difference between CTR and the World Ocean Atlas (WOA), and in the right column the difference between NEWS and CTR is shown. All results are shown in mmol m⁻³.

3.1.2 Denitrification and nitrogen fixation

Denitrification is known to be the main sink for fixed N in the ocean (Gruber, 2004; Codispoti et al., 2001). It occurs both in marine sediments and in the water columns under suboxic conditions, for example in the simulated Bay of Bengal. As a result of these dynamics, if N is added via river discharge, UVic simulates globally higher water column and benthic denitrification rates (Table 3). Note that in these simulations global benthic denitrification and global water column denitrification amount to similar magnitudes, indicating the importance of both in the global N cycle (Somes et al., 2016, 2013). For both processes estimates vary considerably: for water column denitrification estimates are between 50 and 150 Tg N yr⁻¹; for benthic denitrification they are between 100 and 300 Tg N yr⁻¹ (Galloway et al., 2004; Gruber, 2004; Bohlen et al., 2012; Somes et al., 2013). In all our simulations, denitrification rates stay in the range assumed for a balanced fixed-N budget in the pre-industrial ocean (e.g., Somes et al., 2013). However, benthic denitrification is more evenly distributed than water column denitrification. In order to study regional effects it is helpful to include both processes. Nevertheless, models can have a balanced nitrogen budget without including benthic denitrification.

While the global pattern of denitrification is very similar in the simulations with additional riverine N compared to CTR, in proximity to river discharge points, total denitrification rates are higher by up to 1 mol N m⁻² yr⁻¹ (Fig. 9). Somewhat off the coasts however, total denitrification appears lower in the simulations with riverine nutrient supply (by up to 50 mmol N m⁻² yr⁻¹). At the same time, total global N₂ fixation rates decrease in all three simulations compared to CTR (Table 3). Nitrogen fixation is a significant process in the marine nitrogen cycle and a major source of nitrogen in the open ocean. Nitrogen-fixing organisms are able to convert dissolved nitrogen gas (N₂) into ammonia but are limited in their growth by phosphate and iron (Deutsch et al., 2007; Moore and Doney, 2007; Karl et al., 1997; Redfield et al., 1963), and phosphate is not altered by the additional N from river supply. The global rate and geographical distribution of nitrogen fixation are still uncertain. Observations remain sparse and highly variable in space and time. Combined with insufficient understanding of the controls of marine N₂ fixation, this results in high uncertainties in the global pattern of marine nitrogen fixation (Wang et al., 2019; Landolfi et al., 2018; Somes et al., 2013). Deutsch et al. (2007) and Luo et al. (2012) estimated a global nitrogen fixation rate of 140 Tg N yr⁻¹, and most recent studies stay in this range, although some studies suggest that the global rates could be much higher (Wang et al., 2019; Landolfi et al., 2018; Somes et al., 2013; Karl et al., 2002).

The global rates calculated from our experiments with UVic (Table 3) are also higher than the estimates from Deutsch et al. (2007) and Luo et al. (2012). Although previous studies with UVic have given rates of N₂ fixation between 128 and 150 Tg N yr⁻¹ (Landolfi et al., 2017; Keller et al., 2012), the CTR simulation in our configuration estimates global N₂ fixation rates of 219 Tg N yr⁻¹. In our case, this is due to the additional integration of benthic denitrification, which has not always been considered in previous UVic studies. The additional N sink in the form of benthic denitrification promotes conditions that favor N fixers, i.e., diazotrophs, leading to higher nitrogen fixation rates.

In the UVic CTR simulation, N₂ fixation is mostly confined to the tropical and subtropical oceans and is especially concentrated in the northern Indian Ocean, the eastern Pacific Ocean, and the eastern Atlantic Ocean (Fig. 8a). This is comparable to the distribution in Keller et al. (2012) and Somes et al. (2010a), both using UVic in different configurations. The patterns of N₂ fixation are therefore consistent with observations, as far as they are known, with the same limitations as for Keller et al. (2012) and Somes et al. (2010a). For example, in the subtropical North Atlantic, where some of the highest rates of N₂ fixation have been measured (Capone et al., 2005), UVic simulates almost no N₂ fixation at all. The simulation NEWS, DIN+DON, and 2×DIN show that adding riverine N leads to a net decrease in N₂ fixation in nearly the whole area where it occurs but especially near the river mouths (Fig. 8b–d). The main regions where N₂ fixation is significantly decreasing are the Gulf of Guinea, the Bay of Bengal, and near the Amazon River mouth.

In a previous study with UVic, Somes et al. (2016) have shown that increasing atmospheric N deposition could lead to a reduction in N₂ fixation, due to non-nitrogen-fixing phytoplankton being more competitive than N fixers when key nutrients like iron and phosphate are limiting. Here, it is the input of riverine nitrogen that stimulates the reduction in N₂ fixation locally, where N reaches the ocean. Reductions in N₂ fixation can then partly explain the lower NO₃ concentrations at the surface of the tropical and subtropical oceans in NEWS, even though these areas are far from riverine N input (see Sect. 3.1.1 and Fig. 6, first row). Note that these results show the distribution at a steady state after 10 000 years of riverine nitrogen supply. Not all fixed nitrogen is consumed by biological activity, but part of the additional N is also transported with ocean circulation and can “replace” N from nitrogen fixation in regions far off the coast, leading to decreasing N₂ fixation at the surface of the tropical and subtropical oceans.

Figure 8 (a) Global distribution of N₂ fixation in CTR in mmol N m⁻² yr⁻¹. The white areas show regions where rates are smaller than 1 mmol N m⁻² yr⁻¹. (b) Difference in annual vertically integrated rates of N₂ fixation calculated from UVic NEWS and CTR in mmol N m⁻² yr⁻¹. (c) Difference in N₂ fixation in the simulation with riverine DIN and DON. (d) Difference in N₂ fixation in the simulation with twice the amount of DIN. The white areas show regions where differences are smaller than 0.01 mmol N m⁻² yr⁻¹. Local minima can be found near the Amazon River basin (from −356 mmol N m⁻² yr⁻¹ in NEWS − CTR to −382 mmol N m⁻² yr⁻¹ in DIN − CTR), in the Bay of Bengal (from −347 to −646 mmol N m⁻² yr⁻¹ in 2×DIN − CTR), and in the Gulf of Guinea (from −180 mmol N m⁻² yr⁻¹ in NEWS − CTR to −303 mmol N m⁻² yr⁻¹ in 2×DIN − CTR).

Figure 9 (a) Global distribution of total denitrification in CTR in mmol N m⁻² yr⁻¹. The white areas show regions where rates are smaller than 100 mmol N m⁻² yr⁻¹. (b) Difference in annual vertically integrated rates of total denitrification calculated from UVic NEWS and CTR in mmol N m⁻² yr⁻¹. (c) Difference in total denitrification in the simulation with riverine DIN and DON. (d) Difference in total denitrification in the simulation with twice the amount of DIN.

3.1.3 The N-cycle feedback mechanisms

The interaction between N₂ fixation, denitrification, and riverine nitrogen supply can also explain the significant loss in NO₃ in some regions identified before: the Gulf of Guinea, the Bay of Bengal, and the western coast of Central America (Figs. 5 and 7). In addition, these three regions also have in common that they are known to have very low oxygen concentrations. Note that the global volume of ocean minimum zones, defined here as regions with O₂ concentration lower than 70 µmol kg⁻¹, is increasing with higher nitrogen supply, from 52×10⁶ km³ in CTR to 54×10⁶ km³ in NEWS, 56×10⁶ km³ in DIN+DON, and 58×10⁶ km³ in 2×DIN. In the Bay of Bengal, oxygen concentrations, even though higher at the surface in NEWS than in CTR, are very low in the NEWS simulations in the subsurface waters and the whole deeper basin (Fig. 7). These suboxic waters are furthermore located in proximity to riverine N input and high denitrification rates (Fig. 10). While total denitrification rates (benthic and water column denitrification) are already quite high in CTR, they are further increased in NEWS, DIN+DON, and 2×DIN in the northern Bay of Bengal, adjacent to the river delta.

Landolfi et al. (2013) found that the negative feedback mechanism between N₂ fixation and denitrification, generally stabilizing the marine N inventory, can turn into a destabilizing positive feedback, generating runaway N loss, if a close spatial association of N₂ fixation and denitrification occurs. This is due to the stoichiometric imbalance created by the combination of these processes. Denitrification occurs in anoxic or suboxic environments, where nitrate or nitrite can be used as a substitute terminal electron acceptor instead of oxygen. Denitrification consumes 7 mol of NO₃ for every mole of organic N provided by N₂ fixation and remineralized anaerobically via denitrification. If more than one-seventh of the organic N provided via N₂ fixation is denitrified, this leads to a net loss of N by more N lost during denitrification than added via N₂ fixation, called a “vicious cycle” by Landolfi et al. (2013).

In the Bay of Bengal, oxygen concentrations appear higher at the surface in NEWS than in CTR by around 1.5 mmol m⁻³ at least in the southern part of the bay, which could be due to enhanced production. However, compared to the WOA, the oxygen concentrations are still very low in the upper 800 m in NEWS like in CTR, and they are particularly low in the NEWS simulations in the subsurface waters and the whole deeper basin (Fig. 7). The vicious cycle is triggered here by the input of new N from riverine export near oxygen minimum zones, explaining the NO₃ deficit found in the simulated Bay of Bengal (Fig. 5). Note that UVic, similarly to most other biogeochemical ocean models, misplaces the main oxygen minimum zone from the Arabian Sea to the Bay of Bengal (Séférian et al., 2020). In reality, high water column denitrification has been observed in the Arabian Sea, while in the Bay of Bengal highly variable oxygen concentrations seem to inhibit denitrification (Johnson et al., 2019; Bange et al., 2005).

At the end of the simulation, the global marine N inventory is higher by 5278 Tg N in NEWS compared to CTR, which corresponds to 1.12 % of the global N inventory in CTR and 2.3 % of the total riverine N input over the 10 000 years of the simulation. Even for the highest scenario (2×DIN), the total increase in global N represents only +2.53 % of the reference N inventory. Most of the additional N input through river discharge is thus compensated for by the feedbacks of the N cycle. However, relative to the total additional input, the N increase in 2×DIN is higher than in NEWS (+2.6 % compared to +2.3 %), which means that the negative feedbacks are not compensated for in 2×DIN as much as in NEWS. A possible reason for this result could be that the main positive N-loss feedback becomes smaller or that the positive feedbacks, resulting in loss of N, take place in very localized low-oxygen areas that can not expand further (e.g., Bay of Bengal, Amazon River mouth), while riverine N is supplied through river mouths scattered over the world.

Figure 10 Total (water column and benthic) denitrification in the Bay of Bengal. (a) Total denitrification from CTR in nmol N m⁻² s⁻¹. In the white areas total denitrification rates are smaller than 1 nmol N m⁻² s⁻¹. (b) Difference in denitrification between NEWS and CTR. (c) Difference in denitrification between DIN+DON and CTR. (d) Difference in denitrification between 2×DIN and CTR.

Figure 11 N₂ fixation in the Bay of Bengal. (a) N₂ fixation from CTR in nmol N m⁻² s⁻¹. The white color indicates regions where N₂ fixation is smaller than 2 nmol N m⁻² s⁻¹. (b) Difference in N₂ fixation between NEWS and CTR. (c) Difference in N₂ fixation between DIN+DON and CTR. (d) Difference in N₂ fixation between 2×DIN and CTR.

Table 3 Global nitrogen sources (river supply, N₂ fixation) and sinks (denitrification) averaged over the last 100 years of the simulations. All fluxes are given in Tg N yr⁻¹. Note that the global sums from sources and sinks do not exactly add up to zero due to natural variability in the modeled N cycle.

3.2 Marine primary production

The rates of simulated annual global net primary production (NPP) compare well to present-day estimates of annual global NPP (43.5 to 67 Pg C yr⁻¹) derived from satellite measurements (Buitenhuis et al., 2013; Westberry et al., 2008; Carr et al., 2006; Behrenfeld et al., 2005) and vary little between the simulations (Table 4). That is NPP increases only slightly with riverine N supply.

Annually averaged and vertically integrated primary production rates from CTR shows high rates in the equatorial eastern Pacific Ocean, Atlantic Ocean, and Indian Ocean as well as in the upwelling region of the western South Atlantic Ocean (Fig. 12a). This global pattern persists in the simulations with riverine N supply (differences are small in Fig. 12b–d). Nevertheless, with rivers supplying N to the ocean, differences are visible at coastal scales: NPP increases locally, close to the river mouths, especially near the coasts of Europe, China, parts of North America, and parts of South America. Most changes range between −60 and +100 g C m⁻² yr⁻¹. In isolated regions, marine primary production rates are lower in the three NEWS simulations than in CTR. This is particularly striking in parts of the Bay of Bengal but also in the equatorial west Atlantic Ocean to the north of the Amazon River mouth. In large regions of the subtropical and tropical oceans where surface NO₃ concentrations in our simulations are lower than in CTR, NPP differences also show a decrease.

Figure 12 (a) Annual vertically integrated rates of primary production (NPP) calculated from UVic CTR in g C m⁻² yr⁻¹. (b–d) Differences in NPP distribution calculated from different UVic simulations with riverine N export and from UVic simulation without riverine nutrient input (CTR) in g C m⁻² yr⁻¹. (b) NPP from NEWS with riverine DIN. (c) NPP in the simulation with riverine DIN and DON. (d) NPP in the simulation with twice the amount of DIN. The extreme values shown in cyan and magenta in these panels are listed in Table 5.

Behrenfeld and Falkowski (1997)Behrenfeld et al. (2005)Carr et al. (2006)Westberry et al. (2008)Buitenhuis et al. (2013)

Table 4 Mean annual depth-integrated NPP from model data and observations.

Table 5 Minimum and maximum values of NPP difference from CTR simulation as a supplement to Fig. 12.

In NEWS, where only DIN is supplied by the rivers, NPP increases mainly in shelf and near-coastal oceans. The main differences are in the magnitude rather than in the patterns of NPP; however, higher amounts of N added by the rivers in DIN+DON and 2×DIN also impact marine productivity in the open ocean, far away from the river mouths. In the western subtropical and tropical waters NPP decreases with higher N input. These regions correspond to the regions where diazotrophs can be found and where N₂ fixation also decreases in our simulations. In the higher northern latitudes primary production is enhanced off the coastal oceans and in the North Atlantic and the western North Pacific oceans as well as on the Arctic Ocean shelf (Fig. 12d). Two “physical” explanations suggest themselves: first, near the coast the riverine N is consumed until phosphate becomes limiting. Then, the excess N is exported from the coastal oceans, leading to higher productivity farther away. Second, decreasing NPP in the open ocean can be the consequence of a seesaw effect, also called “nutrient robbing”. Because higher N concentrations increase NPP in the coastal ocean, other nutrients, like P, are also consumed here instead of being exported to the open ocean. This can lead to lower rates of primary production farther away (Giraud et al., 2008).

Globally the differences in NPP in NEWS compared to the control simulation are close to the spatial variance of annually averaged NPP in the model (∼ 57 g C m⁻² yr⁻¹). Most changes are smaller than ± 40 g C m⁻² yr⁻¹ (± 2 %), even though locally the changes can be high (down to −122 g C m⁻² yr⁻¹ in the Bay of Bengal or even higher than +500 g C m⁻² yr⁻¹ in the East China Sea). The total increase in NPP varies between +0.7 % (NEWS) and +1.3 % (2×DIN) compared to CTR. These changes reflect in wide parts changes in NO₃ due to the riverine inputs, except for at the higher latitudes and in other regions where light, temperature, or iron limitation occur. In the higher-DIN experiment, NPP is globally a little higher than in the other simulations including CTR, but the distribution shows NPP hot spots near the river mouths, which are compensated for by losses in the subtropical and tropical oceans. In the Bay of Bengal in particular, NPP is enhanced near the outflow of the Meghna River (up to +36 Tg C yr⁻¹), but due to important loss of N further south, total NPP in the basin is lower in all simulations compared to CTR (Table 6).

Table 6 Total NPP, N₂ fixation, and total denitrification from all simulations compared to CTR in some regions with significant decrease or increase in NPP as seen in Fig. 12. NPP is shown in g C yr⁻¹ and N₂ and denitrification in Tg N yr⁻¹.

As stated before, decreasing N₂ fixation together with higher rates of denitrification partly compensates for the additional N from riverine export and thereby buffers the increase in NPP. Note in this regard that a small fraction of NPP is primary production by diazotrophs, which is lower in the NEWS simulations than in CTR (−0.09 Pg C yr⁻¹ for NEWS, not shown here).

In comparison with the NPP pattern in observation (e.g., Behrenfeld and Falkowski, 1997), simulated primary production in UVic is higher in the tropical upwelling regions along the Equator and the northern Indian Ocean, while there is too little productivity on the shelves (Keller et al., 2012). In the simulations with riverine N supply, NPP is increased in the coastal regions and, at least for the 2×DIN experiment, lower in the tropical upwelling regions than in CTR.

The increase in NPP is mainly driven by higher rates near the river mouths, whereas primary production declines in regions where rates of N₂ fixation are lower as a reaction to the input of riverine N, like in the Bay of Bengal and near the Amazon River mouth but also in parts of the open subtropical and tropical oceans (Figs. 8, 11, and 12). This also shows that the response of ocean biogeochemistry depends on the region where riverine DIN first reaches the ocean.

3.3 Simulations with regionally activated riverine nitrogen supply

To answer the question of which rivers have the highest influence on global marine biogeochemistry, we performed five additional experiments with the same configuration as in the NEWS simulation before, but this time with only the rivers of one of five parts of the world transport NO₃ to the ocean. The five scenarios simulate the nutrient supply from North American rivers (NAM), South American rivers (SAM), European and Russian rivers (EUR), Asian rivers (ASIA), and African rivers (AFR).

Figure 13 Difference in the zonal mean concentrations of NO₃ in the main ocean basins (Atlantic Ocean, Pacific Ocean, and Indian Ocean) between regional NEWS and the control simulation (CTR). The regional simulations show export from rivers in North America (a–c), South America (d–f), Europe and Russia (g–i), Asia (j–l), and Africa (m–o). The violet line indicates the 0.0 mmol N contour line. Maximum differences can be found in the Indian Ocean ASIA, with minimum values at −0.9 mmol N m⁻³ and maximum values at 2.5 mmol N m⁻³.

Total riverine N input varies depending on the rivers. Therefore, the amount of N added to the ocean is different in each of the five scenarios. The highest amount of N is added by Asian rivers (11.7 Tg N yr⁻¹). Rivers from South America export 3.5 Tg N yr⁻¹, followed by rivers in EUR (3.2 Tg N yr⁻¹). The lowest-input scenarios are NAM (2.6 Tg N yr⁻¹) and AFR (2.3 Tg N yr⁻¹). The three ocean basins are then affected differently depending on the scenario.

Compared to the control simulation, the differences in NO₃ concentrations are small (0–2 mmol m³), but the patterns differ depending on the origin of NO₃ (Fig. 13). It appears that all ocean basins are most affected by rivers in Europe and Russia (EUR) and least affected by rivers from South America (SAM), although rivers from this region provide the second-highest N supply to the global oceans. The Atlantic Ocean is most affected by rivers in EUR and to a lesser extent in NAM. Asian rivers lead to a local increase in NO₃ concentrations in the North Pacific in the upper 2000 m. Here, NO₃ is trapped because North Pacific waters upwell in the North Pacific. Extra nitrogen is used by local biota and recycled within the North Pacific. But globally and over all depths, rivers in EUR and to a lesser extent in NAM have the biggest impact on NO₃ concentrations in the Pacific. Rivers in SAM slightly decrease NO₃ concentrations in large parts of the Pacific Ocean.

In the Indian Ocean basin, NO₃ concentrations are higher in the simulations NAM, EUR, and AFR. This is because global circulation transports N to remote ocean basins but is not consumed by local primary production, nor does it trigger the vicious cycle described before. In contrast, zonally averaged NO₃ concentrations are lower in the northern Indian Ocean if Asian rivers supply N and thereby enhance denitrification (Fig. 13). N from other regions does not trigger the vicious cycle (Landolfi et al., 2013) in the Bay of Bengal, in the model because it arrives and is used in biological production in other regions before it can reach the Bay of Bengal.

In all simulations with regionally activated riverine nutrient input, N₂ fixation rates are lower than in CTR, with differences ranging between −0.8 Tg N yr⁻¹ in EUR to −7.7 Tg N yr⁻¹ in ASIA (not shown here). This decrease is most prominent in the Bay of Bengal for experiment ASIA because high DIN export from the Ganges Delta gives the advantage to non-N-fixing phytoplankton, which outcompete diazotrophs.

Generally, rivers that enter well-oxygenated eutrophic oceans with little N₂ fixation have the largest impact on the global ocean N inventory. This is especially the case for rivers from EUR and NAM, entering the Atlantic and Arctic oceans at higher northern latitudes. In contrast, the Amazon in SAM is located in an oxygen-deficient region in the tropical Atlantic. The main riverine N supply in ASIA increases N concentrations in the higher northern latitudes of the Pacific but leads to a net loss of N in the Bay of Bengal.

This has consequences for marine productivity: although NPP is higher in NEWS in most of the coastal oceans, where rivers export DIN, NPP is considerably lower in three regions, where the positive vicious-cycle feedbacks dominate – in AFR in the Gulf of Guinea, in ASIA in the Bay of Bengal, and in SAM where the river plume of the Amazon River enters the tropical Atlantic Ocean (Fig. 14).

The global NPP rates range from 54.96 Pg C yr⁻¹ in SAM to 55.08 Pg C yr⁻¹ in the EUR simulation (Table 7). Related to the total amount of N added to the respective simulations via river supply however, EUR contributes considerably more to a widespread increase in NPP, while SAM contributes the least (Table 7). ASIA only increases NPP in a confined part of the Pacific Ocean and mainly decreases primary production in other regions like the Bay of Bengal (Fig. 14).

Figure 14 Difference in the global vertically integrated rates of primary production (NPP) between regional NEWS and the control simulation (CTR). The regional simulations show export from rivers in Africa (a), Asia (b), Europe and Russia (c), North America (d), and South America (e).

According to our simulations, the ecosystem's response to riverine N supply differs depending on the region of the supply and does not always result in an increase in marine NPP. Note that the sum of the NPP changes in the regional experiments is equal to the change in NEWS, except for parts of the Southern Ocean and the eastern Mediterranean Sea. In the Mediterranean Sea, the sum of the regional NPP changes compared to CTR is higher by 1.8 % than the change in global NEWS. This is due to phosphate limitation in this region in the model, limiting increases in NPP in the NEWS simulation.

4 Limitations and further discussion

The results of the simulations with UVic and riverine N have to be evaluated in the context of previous studies. Lacroix et al. (2020) for example found that adding riverine nutrient supply increased primary production essentially in some “hot spots” near the river mouths. While we have observed a qualitatively similar phenomenon in our simulations, these hot spots differ between both studies. This is due partly to the coarser resolution of UVic. The grid configuration used by Lacroix et al. (2020) (GR15) consists of a bipolar grid which resolves the ocean horizontally at around 1.5^∘ and through 40 unevenly spaced vertical layers. Although riverine N is simulated as a fixed percentage of P, dynamic nitrogen fixation by cyanobacteria is included as well as nitrogen deposition and denitrification. From the semi-enclosed seas, which present higher NPP in the study by Lacroix et al. (2020), only the Yellow Sea is adequately resolved in UVic and also shows a significant increase in NPP. But the patterns of primary production differ in several other aspects. Where NPP is substantially increased in the subtropical and tropical eastern Pacific in Lacroix et al. (2020), there is hardly any change in UVic NEWS. In the Bay of Bengal, where we found decreased NPP upon addition of riverine N supply in UVic, their model simulates an increase. The main reason for these differences is the fact that Lacroix et al. (2020) included more than just N from river discharge. Especially the supply in additional phosphate plays an important role in NPP and is for example particularly high in the Bay of Bengal in their simulation. Furthermore, the magnitudes of oceanic nutrient inputs do not differ substantially between the two simulations analyzed in the study of Lacroix et al. (2020): the total N input is 25.2 Tg N yr⁻¹ for the reference simulation (REF) and 27.0 Tg N yr⁻¹ for the simulation with riverine nutrient supply (RIV). The reference distributions of NPP (REF) in Lacroix et al. (2020) also differ with regard to UVic. NPP from UVic is notably higher in the Indian Ocean and the western tropical Atlantic but lower in the Southern Ocean. These higher absolute values of NPP in the open oceans can be accounted for by the parameterization of NPP in our model, where open oceans have to compensate for the lack of higher coastal production in order to achieve estimated annual NPP within the range of 51 to 67 Pg C yr⁻¹ that includes coastal NPP (Keller et al., 2012). But as both reference distributions of NPP differ (from Lacroix et al., 2020, and CTR), it is no surprise that NPP also presents different patterns in the simulations with riverine nutrient supply (RIV and NEWS).

Table 7 Global NPP from different UVic simulations.

Riverine nutrients only reach the ocean in very localized areas. In our simulation with the NEWS data set from Mayorga et al. (2010), we overestimate the effects of adding N from river discharge because DIN is exported directly from the river mouth to the ocean, as our global model does not fully resolve shelf seas and coastal oceans. In reality, a part of these nutrients stays on the shelf or is buried or denitrified in coastal sediments. We also do not account for the buffer effect of the coasts that could be parameterized, as shown by Sharples et al. (2017) and Izett and Fennel (2018). Nevertheless, even without taking these trapping processes into account, the biogeochemical feedbacks of the ocean buffer higher increases in N concentrations. The absolute increase in marine primary production is small (between +0.7 % in NEWS and +1.3 % in 2×DIN). However, relative to the amount of N added to the global ocean, primary production increases yearly by 17.5 Tg C per additional 1 Tg N in NEWS (16.0 Tg C per 1 Tg N in 2×DIN). As we have shown, primary production increases mainly near the river mouths where high nutrient loads are injected in shallow ocean areas, creating production hot spots, while only small changes in production have been found in the open ocean.

Other studies with additional N supply also found only a moderate increase in global primary production rates. Da Cunha et al. (2007) for example predict increases in NPP of up to +5 % for the global ocean but using a high-DIN scenario which includes 7.1 Tmol N yr⁻¹ (corresponding to ∼ 100 Tg N yr⁻¹). Da Cunha et al. (2007) also include silicate, iron, and dissolved inorganic carbon but conclude that riverine N may have the higher impact on primary production.

Like Somes et al. (2016), who also simulated a very small increase in NPP upon the addition of N deposition in their model, we found that decreasing N₂ fixation and increasing denitrification act globally as negative feedbacks. Compared to the total amount of N added by the rivers at the end of the simulation, only 2.3 % (NEWS) to 2.6 % (2×DIN) is retained in the global inventory. The feedbacks compensate for much of the nitrogen addition. In regions of low oxygen concentrations, these feedbacks even overcompensate for the external perturbation in terms of riverine N supply, by forming a vicious cycle (Landolfi et al., 2013), consuming more N than provided by the rivers. This is especially the case in the Bay of Bengal. However, we are aware of the fact that UVic, like several other models, currently misplaces the oxygen minimum zone from the Arabian Sea to the Bay of Bengal. It is likely that N supply by Asian rivers would lead to a somewhat larger increase in the oceanic N inventory if the high nutrient input from the Ganges Delta were not to meet the high denitrification zone in the Indian Ocean (Johnson et al., 2019).

Including other nutrients in addition to N, especially P, could change the setting, especially in regions that are phosphate limited. While this is the logical continuation of the current study, the scope of this project was to explore the consequences of locally high N injections on the N cycle and its feedbacks. Furthermore, as rivers supply relatively more N than P to the global ocean, excess N would still be supplied to the coastal oceans (Turner et al., 2003).

5 Conclusions

In this study a new model component was added to the global Earth system model UVic, simulating DIN supply from river discharge. At the end of the 10 000 years of simulations, the N budget reached a new steady state. The main conclusions regarding the N cycle and marine productivity include the following:

• Riverine N added to the coastal ocean is taken up by near-coastal biology but also exported in the deeper ocean and circulated worldwide.

• Despite the continuous addition of N to the system, global marine N concentrations and marine productivity do not increase substantially (+1.12 % N and +0.72 % for NPP) in our simulations.

• In the coastal regions and especially in some hot-spot regions near the river mouths, riverine nitrogen input leads to higher primary production. Globally, NPP rates increase to up to 17.5 Tg C yr⁻¹ per 1 Tg N yr⁻¹ added to the ocean.

• The globally negative feedbacks of the N cycle buffer most of the increase in NO₃ concentrations and in NPP. N₂ fixation decreases promptly after the beginning of the simulations, partly compensating for the additional N at the surface of the ocean, like the N deposition experiments by Somes et al. (2016) and Landolfi et al. (2017). Water column and benthic denitrification are higher compared to a control simulation without riverine N input and play an important role in low-oxygen regions that moreover tend to expand upon the addition of riverine N supply and generate a net N loss.

• In our regional simulations we have shown that NPP can even decrease locally depending on the region where N reaches the ocean. While N from river discharges from North America and Europe (and Russia) is also circulated and exported to the deeper ocean, N from Asian rivers is trapped in the western Pacific or even partly lost via denitrification in oxygen-deficient regions, like it is in the case of the modeled Bay of Bengal.

• The biogeochemical feedbacks of the ocean buffer further increases in global N concentrations and global NPP. Hence, the result also suggests that ocean fertilization with nitrogen alone (as proposed for example by Harrison, 2017) may not have the desired effect. Indeed, simulated carbon export, evaluated at the 122 m level and including all the shelf regions, increases globally by only 0.06 Pg C yr⁻¹ in our NEWS simulation, representing less than 10 % of the amount estimated by Harrison (2017) to be the upper limit of sequestered carbon in the ocean from ongoing fertilization with nitrogen. On short timescales targeted N fertilization might have some impact, but once the vicious cycle has a chance to start, then the findings would probably be the same, as shown, for example, by Somes et al. (2016). But further research would need to be carried out for targeted spatial and temporal N additions at different levels in N-limited regions.

We have found that like atmospheric deposition, river supply of nitrogen is not only relevant for the coastal system but also for marine biology in the global ocean. But while atmospheric deposition provides only N, rivers also supply P to the ocean. Adding P in addition to N in the coastal oceans may change the response of the ecosystem, especially if N limitation is overcome. Tyrrell (1999) stated that nitrate is the “proximate limiting” nutrient in surface waters, the most limiting nutrient for instantaneous growth. Including phosphate, as the “ultimately limiting nutrient”, could change our story on longer timescales. Future research will therefore include model experiments with the combination of riverine nitrogen and phosphorus.