An explicit estimate of the atmospheric nutrient impact on global oceanic productivity

State-of-the-art global nutrient deposition fields are here coupled to the biogeochemistry model PISCES to investigate their effect on ocean biogeochemistry in the context of atmospheric forcings for preindustrial, present, and future 10 periods. PISCES, as part of the EC-Earth model suite, runs in offline mode using prescribed dynamical fields as simulated by the ocean model NEMO. Present-day atmospheric deposition fluxes of inorganic N, Fe, and P into the global ocean account for ~40 Tg-N yr-1, ~0.28 Tg-Fe yr-1, and ~0.10 Tg-P yr-1. Preindustrial atmospheric nutrient deposition fluxes are lower compared to present-day (~51%, ~36%, and ~40% for N, Fe, and P, respectively). However, the overall impact on global productivity is low (~3%) since a large part of marine productivity is driven by nutrients recycled in the upper ocean layer or 15 other local factors. Prominent changes are, nevertheless, found for regional productivity. Up to 20% reductions occur in oligotrophic regions such as the subtropical gyres in the northern hemisphere under preindustrial conditions. In the subpolar Pacific, reduced preindustrial Fe fluxes lead to a substantial decline of siliceous diatom production and subsequent accumulation of Si, P, and N, in the subpolar gyre. Transport of these nutrient-enriched waters leads to strongly elevated production of calcareous nanophytoplankton further south and southeast, where iron no longer limits productivity. The North 20 Pacific is found the most sensitive to variations in depositional fluxes, mainly because the water exchange with nutrient-rich polar waters is hampered by land bridges. By contrast, large amounts of unutilized nutrients are advected equatorward in the Southern Ocean and North Atlantic making these regions less sensitive to external nutrient inputs. Despite the lower aerosol N:P ratios with respect to the Redfield ratio during the preindustrial period, the nitrogen fixation decreased in the subtropical gyres mainly due to diminished iron supply. Future changes in air pollutants under the RCP8.5 emission scenario result in a 25 modest decrease of the atmospheric nutrients inputs into the global ocean compared to present-day (~13%, ~14%, and ~20% for N, Fe, and P, respectively), without significantly affecting the projected primary production in the model. Sensitivity simulations further show that the impact of atmospheric organic nutrients on the global oceanic productivity turns out roughly as high as the present-day productivity increase since preindustrial times when only the inorganic nutrients’ supply is considered in the model. On the other hand, variations in atmospheric phosphorus supply have almost no effect on the 30 calculated oceanic productivity.


Introduction
Marine primary production is a critical component of the global carbon cycle and important for sustaining the habitability on Earth, although vulnerable to environmental changes (e.g., Steinacher et al., 2010). For example, an estimated decline of subarctic productivity has been reported to accompany the warming of the last 150 years (Osman et al., 2019). Global warming induced by greenhouse gas emissions has increased ocean stratification, reducing the supply of nutrients from subsurface 5 waters and inhibiting the growth of phytoplankton in the surface ocean (Behrenfeld et al., 2006). Thus, the role of nutrient supply by atmospheric deposition will likely be more important in a warmer climate. Several studies have documented the importance of primary production on the surface ocean CO2 concentrations (e.g., Falkowski et al., 2000;Gruber, 2004;Gruber et al., 2009;Le Quéré et al., 2013;Smith, 2019) via the carbon uptake and sinking of the particulate organic matter to the deeper ocean. However, significant uncertainties remain in the projected production among state-of-the-art model simulations, 10 which can range between 2 and 20% for the coupled model intercomparison project phase 3 (CMIP3) and CMIP5 models (Fu et al., 2016;Steinacher et al., 2010), mainly due to the different responses of phytoplankton production to changes in water temperature and stratification (Gröger et al., 2013;Laufkötter et al., 2016;Steinacher et al., 2010).
During primary production, the growth of the phytoplankton functional types (e.g., diatoms and nanophytoplankton) results in a newly formed particulate organic matter within the euphotic zone. These processes are limited however by light, temperature, 15 and nutrients' availability. Nutrient inputs to the euphotic upper ocean result from internal ocean dynamics, such as upwelling or external sources, i.e., input by rivers and atmospheric deposition. The effects of riverine inputs have been, however, widely investigated and are found mostly restricted along the coasts or in marginal shelf basins (e.g., Behrenfeld et al., 2006;Gröger et al., 2013;Holt et al., 2012). Hence, the atmospheric deposition is the only external supply that can reach distal open ocean regions far remote from land. 20 Human activities have heavily perturbed the atmospheric chemical composition and thus the nutrient deposition inputs to the ocean (e.g., Mahowald et al., 2017), but their impact on marine biogeochemistry, and consequently on the oceanic carbon-and al., 2020; Matsui et al., 2018). However, the aerosols from natural and combustion sources tend to be deposited in different regions of the oceans. For example, the subtropical North Atlantic Ocean and the Arabian Sea receive the majority of Fe originated from natural dust aerosols, in contrast to the Pacific and Southern oceans where the Fe-containing combustion aerosols play a more important role compared to atmospheric dust (Ito et al., 2019b).
Present-day global atmospheric DFe and dissolved P (DP) deposition fluxes into the ocean are calculated in the range 0.2-0.4 5 Tg-Fe yr -1 (Ito et al., 2019a;Myriokefalitakis et al., 2018) and 0.10-0.17 Tg-P yr -1 (Mahowald et al., 2008;, respectively. Myriokefalitakis et al. (2016) also demonstrated that DOP can contribute to DP more than 50% over the equatorial oceanic regions. Compared to the present day, DP and DFe emissions may have increased by roughly 3 and 6 times, respectively, since 1850 (Myriokefalitakis et al., 2015. Wang et al. (2014) further showed that DP emissions increased due to the extended use of biofuels in the energy production sector in developing countries as well as emissions due 10 to extensive deforestation in South America and Southeast Asia. By contrast, a significant increase for DFe emissions before the 1990s due to coal combustion is demonstrated, followed-up by a decline due to the implementation of air-pollution abatements (Wang et al., 2015b).
The primary production is strongly linked to both nitrogen and iron acquisitions by marine biota. However, surface oceanic nutrient concentrations are strongly impacted by atmospheric deposition on a regional scale. Nitrogen atmospheric inputs were 15 shown to have a significant effect on marine productivity, export production, and carbon uptake in Low-Nutrient-Low-Chlorophyll (LNLC) regions. For example, the impact of N and P atmospheric deposition on strong oligotrophic regions, such as the Eastern Mediterranean, may lead to an increase in present-day primary production by up to 35% (Christodoulaki et al., 2013), resulting overall in a total phytoplanktonic biomass increase by 16% since preindustrial times (Christodoulaki et al., 2016). The global marine primary production rates are currently estimated in the order of 44 -67 Pg-C yr -1 , based on 20 biogeochemistry calculations and satellite-based estimates (e.g., Aumont et al., 2015;Behrenfeld et al., 2005;Gröger et al., 2013;Krishnamurthy et al., 2007Krishnamurthy et al., , 2009Krishnamurthy et al., , 2010Steinacher et al., 2010). Krishnamurthy et al. (2009) also suggested that the simultaneous anthropogenic N and Fe deposition can increase oceanic productivity by ~1.3 Pg-C yr -1 in the current era compared to the preindustrial conditions. The impact of soluble iron deposition on the carbon export efficiency, although highly uncertain, indicates overall that the Southern Ocean is more sensitive in the projected fire emissions changes compared to other 25 oceanic regions (Hamilton et al., 2020). Assuming, however, complete assimilation of anthropogenic nitrogen by carbon fixation, a new marine biological production up to 0.3 Pg-C yr -1 could be supported (Duce et al., 2008). Although P deposition may account for only a small fraction of production (Krishnamurthy et al., 2010), an increase in both Fe and P deposition can enhance the N2-fixation in LNLC oceans (Mahowald, 2011;Moore et al., 2013b).
A large portion of the global ocean, especially the subtropical gyres, is depleted in nitrate and phosphate, and consequently, 30 sustain low productivity (e.g., Moore et al., 2013a). About 40% of the global ocean is estimated to be N-limited (Krishnamurthy et al., 2009;Wang et al., 2015a), with most of the remaining to be Fe-limited. The relatively larger increase in N than P deposition in many oceanic regions of the globe causes shifts from N to P limitation (Krishnamurthy et al., 2009;Moore et al., 2013a). On the other hand, many studies suggest that anthropogenic Fe deposition is the most important factor for carbon uptake (Krishnamurthy et al., 2010;Okin et al., 2011), mainly due to its positive effect on productivity in High-Nutrient-Low-Chlorophyll (HNLC) regions. Accordingly, the essential role of iron in oceanic productivity is currently well established  and routinely included in marine biogeochemistry models (Aumont and Bopp, 2006;Hajima et al., 2019;Hamilton et al., 2020;Ito et al., 2019b;Moore et al., 2001;Tagliabue et al., 2014Tagliabue et al., , 2016. Another effect of iron is that it stimulates the nitrogen fixation (Camarero and Catalan, 2012;Schulz et al., 2012), because N2-5 fixing species (diazotrophs) have elevated Fe requirements (Kustka et al., 2002). For example, the N2-fixation is found to be suppressed due to iron-limitations in the eastern tropical Pacific (Wang et al., 2019). Hamilton et al. (2020) further demonstrated that changes in iron deposition fluxes into the global ocean may affect up to 70% the marine nitrogen cycle via increases in denitrification and nitrogen fixation rates. Present-day global nitrogen fixation is currently estimated in the range of ~111-163 Tg-N yr -1 (e.g., Aumont et al., 2015;Krishnamurthy et al., 2009;Wang et al., 2019). Wang et al. (2019) suggested 10 that roughly half of the export production in the subtropical gyres is due to the nitrogen from microbial fixation and external inputs, such as rivers and atmospheric deposition. Regarding the atmospheric inputs, however, the aforementioned study demonstrated based on inversion calculations a potential decrease (~10%) of N2-fixation, as a response to the elevated presentday nitrogen emissions.
The present study aims to analyze the impact of a comprehensive representation of atmospheric inputs to the oceanic 15 productivity. For this, a state-of-the-art marine biogeochemical model is used to integrate the recent knowledge of the atmospheric nutrient deposition fluxes into the ocean, driven by natural and combustion emissions, along with further processing during atmospheric transport. The outlined variable composition and varying sources of the deposited nutrients (i.e., N, Fe, and P) used in this work, have been recently modeled with a state-of-the-art atmospheric chemistry and transport model based on preindustrial, present, and future emissions. The description of the biogeochemical model and the 20 parameterizations used in the atmospheric chemistry transport model, which determines the atmospheric deposition fields of this work, are presented in Sect. 2. A detailed description of the regional changes in deposition fluxes and the linked atmospheric processes controlling them is also provided. In Sect. 3, the modeled nutrient oceanic concentrations are presented, and the relevant biogeochemical processes, such as the nitrogen fixation and the primary production, are discussed and compared to estimates from observations and other modeling studies. The role of present-day air pollutants on nutrients' 25 atmospheric deposition is also discussed, via comparison of experiments forced from atmospheric inputs of preindustrial and projected anthropogenic and biomass burning emission scenarios. The impact of the atmospheric nutrients' organic fraction on the global oceanic productivity is assessed in Sect. 4. Moreover, the implications of our findings concerning the above biogeochemistry parameters are discussed in Sect. 5. Finally, the main conclusions are summarized in Sect. 6.

Model description 30
The state-of-the-art biogeochemistry model PISCES (Aumont et al., 2015), enabled within the framework of the European Community Earth System Model EC-Earth (http://www.ec-earth.org/), is here used in offline modus to investigate the impact of atmospheric deposition fluxes of N, Fe, and P on the marine productivity. PISCES (Pelagic Interactions Scheme for Carbon and Ecosystem Studies v. 2), as a part of the Nucleus for European Modelling of the Ocean (NEMO), includes a detailed representation of the lower trophic levels of marine ecosystems. The model simulates the biogeochemical cycles of carbon and the main nutrients (N, P, Fe, and Si), assuming a constant Redfield ratio (i.e., C:N:P = 122:16:1) in organic matter and living biomass. External nutrient sources from atmospheric deposition, rivers, sea ice, sediment dissolution, and hydrothermal vents 5 are also considered. PISCES includes two types of phytoplankton, namely calcareous nanophytoplankton and siliceous diatoms, and it simulates the chlorophyll concentrations and the phytoplankton growth based on the dissolved nutrients' availability (i.e., DP, DN, and DFe for nanophytoplankton and DP, DN, DFe, and DSi for diatoms), temperature and light.
Phytoplankton can be grazed by zooplankton or enters directly into the detritus pool. All particulate organic matter sinking to the bottom undergoes remineralization and the nutrients formerly incorporated during photosynthesis are released again. Thus, 10 PISCES simulates the full inorganic carbon cycle including the biological and the carbonate counter pump. At the ocean surface, air-sea gas exchanges for carbon dioxide, oxygen, and nitrogen are parameterized following Wanninkhof (1992). The model has been successfully tested against the response of oceanic productivity to dust (e.g., Guieu et al., 2014) and climatic variability (e.g., Schneider et al., 2008).

Physical Ocean forcing
The dynamical physical outputs used to force PISCES for this study were produced by the physical ocean model NEMO, following the OMIP (Ocean Modelling Intercomparison Project; Orr et al., 2017) protocol. OMIP aims at harmonizing forcing fields of boundary conditions, as well as validation and analysis procedures among different ocean models. Atmospheric forcing fields are from the CORE II forcing (Coordinated Ocean-ice Reference Experiments -Phase II; Large & Yeager, 2009). 20 CORE II provides a 62-year interannual forcing for the period 1948-2009. The physical model is initialized with gridded observational data from the World Ocean Atlas 2013 Zweng et al., 2013) and then ran for 310 years by repeating the 62-year CORE II forcing. The necessary physical variables to force the offline biogeochemical model PISCES (see Table S1 in the Sup. Mat.) were taken from the last 62-year iteration. To avoid, however, any long-term trends from the spin-up, the multi-year (1948-2009; i.e., the 5 th iteration of the 310-year run) mean of daily forcing fields was calculated. The 25 resulting mean 1-year forcing, thus, contains the mean seasonal cycle and is applied (repeatedly) to drive all simulations with the biogeochemical PISCES offline model. All biogeochemical simulations are initialized and forced with the same physical fields from the average 1-year forcing derived from the OMIP run. Thus, all PISCES offline simulations are drift-free in physical variables. More details of the OMIP protocol can be found in Orr et al. (2017) and a first validation of the OMIP run is provided by Skyllas et al. (2019). 7

The ocean biogeochemistry model
For this study, PISCES uses a ~1 o horizontal resolution with a latitudinal grid refinement in the tropics and 75 layers for the ocean (i.e., ORCA R1), and a timestep of 2700 sec. For the initialization of the ocean biogeochemical fields, the climatological fields of oxygen, nitrate, silicate, and phosphate from the World Ocean Atlas 2009 (WOA; Garcia et al., 2010aGarcia et al., , 2010b along with dissolved inorganic carbon (DIC) and alkalinity from the Global Ocean Data Analysis Project (GLODAP; Key et al., 5 2004) are adopted. The default PISCES configuration (Aumont et al., 2015) uses yearly resolution N deposition fluxes of ~67 Tg-N yr -1 into the global ocean (Duce et al., 2008) assuming that all deposited N into the ocean to be dissolved. Respectively, the Fe, P, and Si atmospheric inputs were calculated from the same (monthly resolution) dust deposition field. Specifically, for the Fe atmospheric input, the Fe content in dust was set to 3.5%, and its soluble fraction was derived based on the simulated monthly resolved Fe solubility fields (Luo et al., 2008;Mahowald et al., 2009) overall, resulting in a soluble Fe input to the 10 ocean of ~0.15 Tg-Fe yr -1 . For the P atmospheric input, the P content in dust was set globally to 750 ppm and with a constant solubility of 10% (Mahowald et al., 2008), resulting in a DP deposition flux of ~0.02 Tg-P yr -1 in the global ocean. Finally, for the Si atmospheric input, a constant fraction of 30.8% in dust was set, assuming that 7.5% of the deposited total Si as soluble (~6.35 Tg-Si yr -1 ) and, thus, upon deposition entered in the dissolved silicate pool of the model.
In contrast to previous studies (e.g., Aumont et al., 2015), the new N, Fe, and P atmospheric deposition fields considered here 15 (see Sect. 2.2) are all calculated based on: 1) emissions of natural and combustion nutrient containing aerosols, 2) detailed atmospheric gas-and aqueous-phase chemical schemes, and 3) mineral dissolution processes due to atmospheric acidity and organic ligands in aerosol water and cloud droplets.
Note that, as for the default PISCES configuration, the Si deposition fluxes into the ocean are only based on the new dust 20 deposition fields coupled in the model. Moreover, for this work, no extra optimizations for the iron scavenging parameters have been applied, since the default PISCES configuration already considers a variable iron solubility on the dust deposition inputs (Aumont et al., 2015). The simple chemistry scheme available in PISCES is here used, which is based on one ligand (L) dissolved inorganic Fe and one dissolved complexed iron (FeL). The ligand concentration in the ocean is kept constant, equal to 0.6 nmol L −1 and the scavenging rate by dust is equal to 150 d −1 mg −1 L (see Aumont et al., 2015 and ref. therein). For 25 clarity, we note that the ocean and biogeochemistry modules used for this study may slightly differ from the version recently used in EC-Earth CMIP6 simulations since at the time the simulations were carried out the final version of EC-Earth was still not released.
Two transient simulations from 1651 to 2100 are here performed to study the impact of nutrients' atmospheric input on global marine productivity: 30 1) a standard (STD) simulation accounting for the inorganic fractions of the deposited atmospheric nutrients (N, P, and Fe) into the global ocean, and 2) a sensitivity (ORG) simulation, as for STD but also accounting for the organic fractions of the deposited atmospheric nutrients (N, P, and Fe).
We here present results for the preindustrial (PAST: 1851-1870 average), present-day (PRESENT: 2001-2020 average), and future projected (FUTURE: 2081-2100 average) periods. For all PISCES simulations, the first 200 yrs. (i.e., 1651-1850) are not interpreted but considered as a spin-up to reach a quasi-equilibrium state in the model with a well-ventilated upper ocean. 5 Moreover, the atmospheric CO2 mixing ratio is set to the preindustrial value of 284.7 ppm, to effectively isolate the impact of atmospheric deposition on the marine biogeochemistry parameters. To account, however, for potential drifts in the deeper ocean layers a control (CTRL) simulation, as for STD but using only preindustrial (i.e., the year 1850) atmospheric nutrients (N, P, and Fe) inputs into the global ocean, is also performed. Figure S1 (see Sup. Mat.) demonstrates that for the main ocean basins the drift in vertically integrated primary production is minimal and clearly below the signal imposed by the altered 10 nutrient deposition after 1850. This holds even for the Southern Ocean where the impact of atmospheric deposition is typically weak due to the absence of neighboring emission sources. Nevertheless, all model results presented in this work have been adjusted by subtracting the drift of the control run from STD.

The atmospheric nutrient inputs
All atmospheric nutrient inputs coupled to PISCES are derived from the offline global atmospheric Chemistry-Transport 15 Model (CTM) TM4-ECPL. The CTM is driven by the ECMWF (European Centre for Medium-Range Weather Forecasts) Interim reanalysis project (ERA-Interim) meteorology (Dee et al., 2011) for the year 2010, and uses a horizontal resolution of 3 o in longitude by 2 o in latitude, with 34 hybrid layers up to 0.1 hPa. The CTM simulates the gas-phase chemistry along with the major non-methane volatile organic compounds, as well as all major aerosol components, such as dust, sea-salt, organic aerosols (OA), black carbon (BC), SO4 2-, NH4 + , and NO3 -. The thermodynamic equilibrium model ISORROPIA II (Fountoukis 20 and Nenes, 2007) is used to estimate the water content and the acidity of hygroscopic aerosols, accounting also for the impact of crustal materials (i.e., Ca 2+ , Mg 2+ , K + , Na + , Cl -) from mineral dust and sea-salt (Myriokefalitakis et al., 2015). The in-cloud acidity is mainly controlled in the model by strong acids, i.e., sulphuric acid H2SO4/SO4 2-, methanesulphonic acid (MSA/MS -), HNO3/NO3 -, and bases (NH3/NH4 + ), accounting also for the dissociation of hydrated CO2, sulfur dioxide (SO2), and oxalic acid (Myriokefalitakis et al. 2011). The CTM further considers the multiphase chemistry secondary aerosol production in cloud 25 droplets and aerosol water (Myriokefalitakis et al., 2011), as well as the secondary organic aerosol (SOA) formation via gasto-particle partition over land and oceanic regions (e.g., Myriokefalitakis et al., 2010). The anthropogenic (including ships and aircraft emissions) and biomass burning emissions from the historical Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) database (Lamarque et al., 2013) are used in the CTM, and the RCP8. Note that we use here nutrient atmospheric deposition fields based on simulations with the latest updates of the CTM, as recently published by Kanakidou et al. (2020), which are all based on an improved on-line dust emission scheme (Tegen et al., 2002). As a result, the global dust source for all simulations is here equal to ~1287 Tg yr -1 for the year 2010, well comparable to a multimodel estimate of ∼1257 Tg yr −1 as reported by Huneeus et al. (2011). For example, the global annual mean Si deposition input rate to the global ocean, which is solely related to the dust deposition, is equal here to ~4.7 Tg-Si yr -1 and overall is in the range of other estimations when a solubility fraction of 7.5% is considered (e.g., see Krishnamurthy et al., 2010). Nevertheless, some small differences on the total amount of deposited nutrients over the ocean compared to previously published results (e.g., Kanakidou et al., 2016;Myriokefalitakis et al., 2015Myriokefalitakis et al., , 2016 are expected, due to the various updates of 5 the code, the different year of simulation, or even the different definitions of the oceanic regions due to the applied horizontal analysis in PISCES (i.e., the ORCA R1).
Simulations with the atmospheric transport and chemistry model are, nonetheless, extremely expensive. Therefore, limitations in available computational resources made it necessary to reduce the CTM simulations to representative single years for 1) the preindustrial state (before 1850), 2) the present-day state (representing the year 2010), and 3) a mid-century (2050), as well as, 10 an end of century (2100) state. However, as the typical residence time of tropospheric aerosols is in the order of days, the atmospheric depositional fields used in PISCES represent a well equilibrated atmospheric chemistry and deposition flux, without the need of time transient simulations.
For the ocean biogeochemistry model spin-up (i.e., from 1651 to 1850) the preindustrial field (the year 1850) was applied.
After the 200 years spin-up period, the atmospheric deposition input data for the STD and ORG simulations were linearly 15 interpolated from preindustrial to present-day conditions (i.e., the year 2010) to smoothly capture the transition from past to the modern conditions (e.g., Krishnamurthy et al., 2009). Respectively, the deposition data from the present day were linearly interpolated to the projected estimates (i.e., the years 2050 and 2100). Note that for all temporal and spatial interpolations of this work (as well as, for the drift corrections), the Climate Data Operators (CDO v.1.9.8) software, as provided by the Max Planck Institute for Meteorology, is here used (https://code.mpimet.mpg.de/projects/cdo/embedded/cdo.pdf; last access 20 29/02/2020). An example of the globally averaged N, Fe, and P atmospheric deposition data as simulated by the CTM and applied in PISCES is presented in Fig. S2. Overall, the here discussed simulations should be considered as idealized sensitivity experiments to estimate the response on the ocean surface biogeochemical properties to changed atmospheric deposition.

Nitrogen
For the calculation of the atmospheric nitrogen deposition fluxes, the CTM uses primary emissions of NOx, NH3, marine 25 amines, and emissions of particulate organic nitrogen (ON) from natural and anthropogenic sources (Kanakidou et al., 2012(Kanakidou et al., , 2018. The particulate ON is linked to the OA tracers using varying N:C molar ratios (Kanakidou et al., 2012), as well as to the SOA formation under high NOx-to-VOC conditions (Myriokefalitakis et al., 2010;Tsigaridis et al., 2006). Amines of marine origin in the gas phase are also considered to form amine salts (Myriokefalitakis et al., 2010). A more detailed description of the N-cycle parameterization in the CTM can be found in Kanakidou et al. (2016Kanakidou et al. ( , 2018. 30 Figure 1b presents the annual mean spatial distribution of dissolved nitrogen deposition as considered in PISCES for the STD simulation. The present-day DIN (oxidized and reduced dissolved inorganic nitrogen) deposition fluxes into the global ocean are estimated to be ~40 Tg-N yr -1 (Table 1). DIN deposition shows the highest fluxes downwind industrial areas of the Northern Hemisphere and the tropical biomass burning regions due to the enhanced NOx emissions, as well as downwind Europe, China, and Indonesia, reflecting overall the high strength of NH3 emissions in these regions. DIN deposition exceeds 1 . 10 -3 kg-N m -2 yr -1 downwind the eastern United States, Europe, India, China, and Indonesia (Fig. 1b). Some low DIN deposition fluxes over the remote ocean are related to the recycling of the marine NH3 sources. Moreover, when the ON is accounted for the dissolved nitrogen deposition, as considered in the ORG simulation ( Fig. S3a), much higher nitrogen deposition fluxes are calculated in 5 the tropics mainly due to the large contribution by primary biogenic particles, the biomass burning emissions, along with the ON production during SOA formation (Fig. S3b). For comparison, we note that the total river supply of bioaccessible nitrogen in the model is roughly 36 Tg-N yr −1 (Aumont et al., 2015). Note, however, that the total present-day dissolved nitrogen deposition estimate (~58 Tg-N yr -1 for the ORG simulation; Table 1) is lower compared to the N deposition fluxes (~67 Tg-N yr -1 ) as taken from Duce et al. (2008) and used in previous PISCES configurations (Aumont et al., 2015). 10 Compared to the present day, almost all ocean basins (except some parts of the South Indian and the South Pacific Ocean) display a substantially lower (>50%) nitrogen deposition flux during the preindustrial era (Fig. 1a). Inorganic nitrogen inputs to the ocean have significantly increased for present-day along the coasts of the African, Australian, and the South American continents, downwind densely populated areas in the northern hemisphere, such as the east coast of North America and Europe, as well as regions with intensive agriculture downwind the coast of East and South Asia . Globally, 15 present-day atmospheric inorganic nitrogen inputs to the ocean have increased by a factor of ~2 since 1850 (Table 1) due to the respective increase of NH3 and NOx emissions. The projected atmospheric DIN inputs to the ocean (Fig. 1c) indicate also a decrease (~15%), although much less significant compared to preindustrial times (Table 1), since the reduction in NOx emissions is projected to be compensated by the continuing increase in NH3 emissions. Note that the preindustrial ON deposition fluxes into the ocean are estimated to be rough of the same magnitude as the inorganic nitrogen oceanic input (i.e., 20 ~15 vs. ~20 Tg-N yr -1 ; Table 1). Projections under the RCP8.5 scenario, however, indicate that NHx emissions will gain in importance, resulting overall in a weaker contribution of the oxygenated inorganic nitrogen to the atmospheric deposition into the ocean for future conditions.

Iron
The global atmospheric Fe deposition fluxes are parameterized in the CTM considering primary Fe emissions associated with 25 minerals in dust and combustion processes. The Fe content of dust minerals is based on detailed mineralogy maps (Nickovic et al., 2012), accounting also for an initial soluble Fe content in the mineral emissions (Ito and Xu, 2014), overall resulting in a mean Fe content of about 3.2% in dust emissions. For the Fe-containing combustion aerosols, the CTM accounts for emissions from biomass burning, coal, and oil combustion (Ito, 2013;Luo et al., 2008) with dissolved Fe contents of 12%, 8%, and 81%, respectively. The CTM further accounts for acid-and organic ligand-solubilizations of dust aerosols, in both 30 aerosol water and cloud droplets, as well as, for the aging (i.e., the conversion of insoluble to soluble) of the Fe-containing combustion aerosols via atmospheric processing. More details on the atmospheric Fe-cycle set-up can be found in Myriokefalitakis et al. (2015Myriokefalitakis et al. ( , 2018, along with updates from Kanakidou et al. (2020). Figure 1e presents the annual mean spatial distribution of dissolved iron atmospheric deposition fluxes, as considered in PISCES for the STD simulation. DFe deposition fluxes into the ocean present strong spatial variability (Fig. 1e) with an annual flux of ~0.28 Tg-Fe yr −1 for the STD and ∼0.35 Tg-Fe yr −1 for the ORG simulation (Table 1). Both estimates are well in the range of the mean DFe deposition flux into the ocean (0.2-0.4 Tg-Fe yr −1 ), as derived from the GESAMP model intercomparison study (Myriokefalitakis et al., 2018). For comparison, we note that the total riverine Fe supply in the model 5 equals 1.45 Tg-Fe yr −1 (see Aumont et al., 2015 and ref. therein). The highest annual mean DFe deposition fluxes into the ocean occur downwind dust source regions. Significant deposition rates are also found in the outflow of tropical biomass burning regions (i.e., Central Africa and Indonesia), reflecting the importance of combustion processes. Annual mean DFe deposition rates of ~10 -5 kg-Fe m −2 yr −1 are considered for the tropical Atlantic Ocean, as well as, for the Indian Ocean under the influence of the Arabian and Indian peninsulas. For the Southern Ocean, the DFe atmospheric inputs (up to ~10 -6 kg-Fe 10 m −2 yr −1 ) are mainly associated with the Patagonian, the Southern African, and Australian deserts. Figure S3c further presents the annual mean DFe deposition fluxes into the ocean when the organic fraction is also considered. The organic-bound Fe is produced in the CTM during the organic-ligand dust dissolution processes (i.e., as Fe(II/III)-oxalato complexes), and further accounted for a fixed 0.1% fraction in Fe-containing combustion aerosols (Kanakidou et al., 2018;Myriokefalitakis et al., 2018). Note that although this estimate is highly uncertain due to the lack of observational data, and it overall appears to 15 contribute modestly to the global DFe atmospheric input to the ocean (~0.07 Tg-Fe yr -1 ), the deposition of organic-bound Fe aerosols can be potentially important (up to ~40%) in the remote tropical Pacific and the Southern Ocean (Fig. S3d) where the atmospheric Fe concentrations are extremely low and are mainly occurring due to biomass burning and anthropogenic (i.e., ship) combustion emissions (Ito et al., 2019a).
The CTM calculates increases in DFe deposition rates since preindustrial times, as stronger Fe combustion emissions and more 20 efficient dust dissolution rates due to a more acidic environment occur in the modern era; overall, accounting for about 1.5 times higher DFe atmospheric input to the global ocean for present-day (Table 1). On the other hand, the derived DFe in the ocean for the future atmosphere is calculated ∼14% lower than present-day conditions (Table 1), owing to the projected changes in anthropogenic emissions and air-quality. For the preindustrial era, the largest differences in atmospheric DFe deposition fluxes compared to present-day are considered in the Northern Indian, the North Pacific, and the tropical North 25 Atlantic Oceans (Fig. 1d). Lower PAST atmospheric DFe deposition fluxes compared to present-day are simulated in remote oceanic regions away from dust plumes, in the tropical and subtropical Pacific Ocean, also due to the increased present-day atmospheric processing. The combustion source of DFe turns out, however, to be rather important near industrial and biomass burning sources, such as downwind South and East Asia, where dust emissions are lower. For future emission scenario, smaller changes are derived, with a general decrease of the atmospheric DFe input in most parts of the global ocean (Fig. 1f), except 30 for some increases in the Eastern North Pacific Ocean and over the very low atmospheric Fe concentrations regions in the remote Southern Ocean.

Phosphorus
The atmospheric P-cycle is calculated based on emissions of insoluble mineral-P, phosphate, and insoluble and soluble organic P (OP), with the resulted DP deposition fluxes in the CTM being driven by natural (i.e., dust, bioaerosols, sea spray and, volcanic aerosols) and combustion P-containing aerosol emissions. Acid solubilization of dust particles (i.e., conversion from mineral P to phosphate) takes place in both aerosol water and cloud droplets, along with the aging of OP-containing aerosols 5 during atmospheric processing. The CTM accounts for two P-containing insoluble minerals (the fluoroapatite and the hydroxyapatite) in dust based on soil mineralogy maps (Nickovic et al., 2012), as well as for OP present in soil's organic matter (Kanakidou et al., 2012). A solubility of 10% is applied to all P-containing dust emissions in the CTM. For P-containing combustion aerosols, the CTM accounts for anthropogenic (i.e., for fossil fuel, coal, waste, and biofuel) and biomass burning emissions, based on observed P/BC mass ratios (Mahowald et al., 2008). Sea-spray and volcanic aerosols account for a rather 10 low DIP global source, in contrast to bioaerosols which are estimated to contribute significantly to DOP. The naturally emitted OP by bioaerosols is overall represented by bacteria, fungi, and pollen (Myriokefalitakis et al., 2017). More details on the Pcycle representation in the CTM can be found in Myriokefalitakis et al. (2016).
The present-day global annual DIP deposition flux into the global ocean for the STD simulation accounts for ~0.10 Tg-P yr −1 (Table 1), presenting also a strong spatial variability (Fig. 1h). The highest DIP deposition input rates occur downwind major 15 dust source regions (~5·10 -6 kg-P m −2 yr −1 ), owing mainly to the phosphate content in dust emissions. High DIP deposition fluxes (~1·10 -6 kg-P m −2 yr −1 ) also occur downwind heavily forested regions (i.e., Amazonia, Central Africa, and Indonesia) due to the enhanced biomass burning sources. The combustion of anthropogenic origin further contributes to the DIP deposition flux into the ocean, such as downwind South and East Asia (~1·10 -6 kg-P m −2 yr −1 ). Notable deposition rates are also illustrated away from dust sources, such as in the Northern Atlantic and Pacific Ocean (~5·10 -7 kg-P m −2 yr −1 ), due to the mineral P 20 solubilization during atmospheric transport and somewhat lower rates occur in the Southern Ocean. Moreover, the consideration of organic fraction to atmospheric DP inputs to the global ocean (Fig. S3e) results in a ~50% increase in the global DP deposition flux (Table 1), but stronger regional increases of up to 80% can be also seen (Fig. S3f). Note also that in PISCES, DP fluxes of roughly 3.7 Tg-P yr −1 are also delivered to the ocean by rivers (Aumont et al., 2015).
An increase in the global DIP deposition of ∼40% is considered here for present-day (i.e., 0.06 Tg-P yr −1 ) compared to 25 preindustrial conditions, and a modest decrease of ∼20% is projected for the year 2100 (i.e., 0.08 Tg-P yr −1 under the RCP8.5 scenario), as shown in Table 1. The regional differences for the preindustrial times appear, however, to be stronger (up to 60 %; Fig. 1g), especially downwind highly populated regions of the Northern Hemisphere, in the Atlantic and Pacific Oceans.
DIP deposition fluxes are also projected to decrease over the midlatitudes of the Northern Hemisphere where human activities dominate under the RCP8.5 scenario, with the largest changes (up to 40%) to be calculated downwind China and Australia. 30 Significant changes (~20 %) are also illustrated in the Arabian Sea and the Bay of Bengal due to the expected increases in population (Fig. 1i).

Oceanic nutrient concentrations
Nitrate: The simulated annual mean nitrate surface concentrations in the seawater for the present-day and the relative differences for past and future eras are presented in Fig. 2. The present-day surface nitrate distribution shows high concentrations along the equatorial divergence where nutrients are upwelled, and solar insolation supports good light 5 conditions throughout the year. In the high latitudes, cooler water temperatures and seasonally damped light conditions reduce nutrient consumption by biological productivity, resulting overall in elevated annual mean nutrient concentrations. High surface concentrations are also calculated in the high latitude Southern Ocean where the deep convection around Antarctica maintains high nutrient transports to the surface and productivity is limited by a thick mixed layer, lower water temperatures, and reduced light conditions. Elevated nutrient concentrations are likewise simulated in the Eastern Equatorial Pacific, and the 10 Subarctic North Pacific, i.e., the well-known HNLC regions. All in all, this reflects a reasonable simulation of the abiotic oceanographic drivers in the model. A comparison between the simulated present-day surface nitrate concentrations and the compiled data from WOA is given in Fig. S4.
Increased atmospheric nitrogen deposition fluxes from the preindustrial to the modern era (Table 1) result in a respective increase in surface nitrate concentrations in almost all oceanic regions, with some exceptions in the Eastern Equatorial and the 15 subpolar Pacific Ocean (Fig. 2a). In remote oceanic areas, far from any coastal or riverine nutrient supply and thus strongly nutrient-limited, the higher present-day inorganic nitrogen and iron atmospheric inputs compared to the preindustrial era (Figs. 1a,d) increase primary production (see also Sect. 3.3). This may, overall, lead to increased export of the surface seawater nitrogen to the deeper ocean in the form of sinking biogenic particles in these oligotrophic oceanic regions (e.g., Krishnamurthy et al., 2007Krishnamurthy et al., , 2010. For the future conditions, however, the model calculates both negatives and positive changes in the surface 20 nitrate concentrations (Fig. 2c), resulting from an overall decrease (~4%) of the global inorganic nitrogen oceanic input (via atmospheric deposition and N2-fixation) (Table 1). Indeed, for the future era, a decrease of atmospheric inputs in almost all oceanic regions is calculated, except for the coastlines of Southeast Asia, Africa, and South America (Fig. 1c). Figure 2 also presents the annual mean surface concentrations of iron for present-day in STD simulation, along with the respective past and future relative differences. Present-day surface iron distribution (Fig. 2e)  For the past (Fig. 2d) and future (Fig. 2f) conditions, the model generally calculates lower iron surface concentrations, 30 reflecting overall the respective decreases in DFe deposition fluxes into the ocean (Figs. 1d,f, respectively). Consequently, the strongest declines are found in the northern hemisphere, especially in the mid to high latitude Pacific and Atlantic. An exception is the NW Pacific where higher iron input in FUTURE (Fig. 1f) results in elevated oceanic iron concentrations. A comparison of the simulated present-day surface oceanic iron concentrations with available DFe oceanic observational data (Tagliabue et al., 2012) is given in Fig. S5.

Phosphate:
The annual mean present-day phosphate surface concentrations as calculated for the STD simulation together with the respective relative differences for past and future eras are presented in Figs. 2g-i. The surface phosphate distribution in the model shows high concentrations along the equatorial divergence where nutrient-rich deep water is upwelled, as well as in the 5 high latitudes, with the highest surface concentrations to be simulated in the Southern Ocean. A secondary maximum is calculated in the Eastern Equatorial Pacific and the Subarctic North Pacific; both regions are subject to large-scale upwelling of deep waters. Note that in general, phosphate concentrations in deep waters are higher than at the surface where nutrients are removed by biotic productivity and exported by sinking particulate organic matter to the deeper ocean. However, due to the constant Redfield Ratio (i.e., C:N:P = 122:16:1) applied in the model (Aumont et al., 2015), the phosphorus cycle is closely 10 related to that of nitrogen, as both are subjected to the same large-scale physical processes and circulation in the ocean. An exception is the N2-fixation that acts as an additional external source for inorganic nitrate in the model. A comparison of the simulated surface phosphate concentrations with the WOA data is given in Fig. S6.
Despite the roughly 1.7 times increase in the phosphate deposition inputs to the ocean from preindustrial to the modern era (Table 1), the preindustrial surface phosphate oceanic concentrations are calculated 20-50% higher in most oceanic regions 15 ( Fig. 2g), except for the Southern Ocean where no significant change is calculated. Accordingly, although the projected decrease of the global phosphate input (Table 1), higher phosphorus surface oceanic concentrations are simulated for the future, up to ~20% (Fig. 2i). The main reason for these elevated phosphate concentrations is the decreased primary production almost everywhere (Figs. 3d,f). Indeed, as nitrogen is the limiting factor for phytoplankton in the open ocean, the primary production rates have been lowered, following the lower nitrogen deposition rates. Accordingly, less phosphate consumption by 20 phytoplankton growth takes place and this outcompetes the effect of lowered P deposition, leading to relatively higher phosphate concentrations. The effect of the decreased productivity on phosphate concentrations would be, however, even stronger if it was not being partly compensated by higher N2-fixation rates. Overall, this points to the marine biogeochemical processes as an essential factor in controlling the phosphorus concentrations at the surface ocean, rather than the depositional fluxes (see also Sect. 3.3.1). 25

Nitrogen fixation
For the STD simulation, the nitrogen fixation is calculated to about 112 Tg-N yr -1 for the present day (Table 1). The respective relative differences compared to past and future periods are also presented in Figs. 3a,c. Compared to modern times, the model calculates a significant decrease (up to 20%) in nitrogen fixation in the tropical and subtropical Pacific and the subtropical Atlantic Ocean for the preindustrial era. On the contrary, downwind Bay of Bengal and Indonesia, nitrogen fixation is higher 30 for the preindustrial conditions, due to the lower nitrogen deposition fluxes accounted for in the model (see Sect. 2.2.1). Finally, the nitrogen fixation rates present very low differences in the Equatorial Pacific, Equatorial Atlantic, and the Southern Indian Ocean (0-10%) away from land sources (Fig. 3a). Note, however, that nitrogen fixation in PISCES is restricted to warm waters (i.e., above 20 o C); therefore, the strong reductions of nitrogen deposition in the mid to high latitude North Pacific in PAST have no direct impact on nitrogen fixation rates. In the subtropical Pacific, the reduced nitrogen fixation rates mainly reflect the diminished iron input (Fig. 1d). On a global scale, the model calculates overall only a small decrease (~0.2%; Table 1) in preindustrial nitrogen fixation rates compared to present-day, mainly as a result of the decreased soluble iron inputs in the subtropical North Pacific (Fig. 1d). For the future conditions, the model likewise calculates a modest decrease in the global 5 nitrogen fixation (~1% ; Table 1) along with decreased iron inputs to the ocean (Fig. 1f), resulting overall in some lower rates of up to 10% in the Equatorial Pacific Ocean (Fig. 3c).

Primary production
The present-day annual mean primary production together with the relative differences to past and future periods are presented in Figs. 3d-f. The primary production distribution in the open ocean shows high rates along the equatorial divergence and in 10 the high latitudes, where nutrients concentrations are high. The decreased nitrogen deposition during preindustrial conditions compared to present-day results in lowered primary production rates almost in all oceanic regions (Fig. 3d). A projected modest decrease of primary production rates is also calculated by the model (Fig. 3f), owing to the lower (~14%) dissolved iron deposition fluxes implemented from present to future conditions (Fig. 1f).
The present-day modeled globally integrated production (~47 Pg-C yr -1 ; Table 1) is, however, lower compared to satellite-15 based estimates from SeaWiFS (Behrenfeld et al., 2005), obtained equal to ~60-67 Pg-C yr -1 , but in the range of estimates from other studies (e.g., 23.9 -49.1 Pg-C yr -1 ; Steinacher et al., 2010). A more detailed comparison between satellite estimates (Behrenfeld et al., 2005) and model simulations of the global primary production is presented in Fig. S7. Overall, the simulated primary production reproduces the main features derived from satellite-based observations (Figs. S7a,d). The model simulates relatively lower rates (~200 mg-C m -2 day -1 ) in the subtropical gyres and higher rates (> 500 mg-C m -2 day -1 ) in upwelling 20 regions, the North Atlantic and the Southern Ocean (Fig. S7g). Nevertheless, as also known from other modeling studies, the primary production in the tropics might be overestimated, whereas in higher latitudes it might be underestimated (e.g., Steinacher et al., 2010). High productivity in the open ocean areas is linked to upwelling areas, such as the equatorial divergence zones or coastal upwellings like at the west coasts of South Africa or northern South America. Primary productivity rates are also underestimated in the North Atlantic region (Fig. S7g), probably due to the low accumulation of nutrients during 25 winter as already discussed. Figure S7 further illustrates the simulated primary production for the boreal winter (DJF) and summer (JJA), compared to the respective observation-based data. The model generally compares reasonably for both seasons; however, the modeled primary production rates (Figs. S7e,f, respectively) are calculated lower in the high latitudes. The pronounced seasonality due to light and temperature limitations of phytoplankton growth in the North Atlantic is well captured by the model. During the warm season, higher water temperatures and better light conditions increase primary production and 30 nutrient consumption by phytoplankton growth in the high latitudes. This results in depleted nutrient concentrations in the North Atlantic during boreal summer. In the North Pacific, which is known as an HNLC region, nutrients remain at higher levels due to insufficient iron support. This is somewhat underestimated in the model as the simulated nitrate concentrations are too low (Fig. S4f), probably related to the slightly higher iron concentrations in the upper 100 m North Pacific (Fig. S5f).
The same biases mentioned for phosphate are likewise seen in the nitrate concentration as most features are mainly controlled by the model physics. Overall, only nitrogen fixation and denitrification can modulate the nitrate concentrations, apart from all other processes that also influence the phosphate in the same way (e.g., productivity, remineralization dissolution, etc.).

Role of depositional nutrient elemental ratios 5
Despite the relatively strong changes in total atmospheric nutrient supply from PAST to FUTURE (Table 1), the impact of atmospheric nutrients on the global productivity rates remains low in the model. This is, nevertheless, not unexpected, as the atmospheric nutrient supply constitutes only a small fraction of the total ocean nutrient inventory. In addition, oceanic regions that are not nutrient-limited today are less sensitive to external nutrient supply. Finally, a large part of primary production is regenerated by remineralized nutrients from particulate organic matter (mainly detritus) in the upper ocean layer. 10 To further identify the oceanic regions that are particularly sensitive to changes in external nutrient inputs, the limiting factors for local productivity in the model are investigated. Figure 4a displays limitations due to nitrogen or phosphorus. High values (indicating low limitation) are seen in regions that are subject to intense upwelling, like in the equatorial divergence zones or the western margins of NW and southern Africa and South America (coastal upwelling). Accordingly, these regions are less sensitive to atmospheric deposition as nutrients are supplied from deep ocean layers. Lower nutrient limitation is likewise seen 15 in the mid to high latitudes where limitations by temperature and light ( Fig. 4b) limit the growth rates. Exceptions are the North Pacific, the Southern Ocean, and the equatorial Pacific where iron limitation matters (Fig. 4c). Consequently, the model's nutrient sensitivity is larger in the subtropics (in particular in the subtropical gyres), where good light conditions and warm waters support high growth rates. Furthermore, these regions are far from land nutrient sources and so, a major part of total primary production relates to regenerated production (i.e., with low rates of external nutrient supply) which is limited by 20 nutrients as both temperature and light are sufficient. This makes productivity in the subtropical gyres sensitive to changes in the external atmospheric nutrient.
In regions with significant macronutrient limitations, the elemental ratio of deposited N:P can be, however, rather important.
To estimate the relative impact of the changes in this ratio, we calculated the modeled nitrogen concentrations relative to the model's Redfield ratio (Fig. 5a). For PRESENT, the model exhibits almost everywhere a deficiency with respect to nitrogen. 25 This is in good agreement with data from the WOA, which likewise indicates a predominant nitrogen deficiency (Fig. 5b).
Next, the N:P ratio relative to the Redfield as supplied by atmospheric deposition for PRESENT together with the changes in PAST and FUTURE is derived (Fig. S8). Overall, a strong excess of N compared to P for modern times is indicated. As a consequence of the model's nitrogen deficiency (Fig. 5a), the atmospheric nitrogen excess maintains higher productivity than without the atmospheric supply. For preindustrial times, however, the atmospheric N:P ratio is almost everywhere reduced, 30 further increasing the N-deficiency. Hence, rather the lowered atmospheric nitrogen inputs than the lowered phosphorus inputs in PAST and FUTURE are responsible for the diminished productivity in these experiments. To further demonstrate this, we carried out an additional sensitivity simulation (namely PIP simulation) where the phosphorus atmospheric deposition fluxes kept constant at preindustrial levels, while the other studied atmospheric inputs (i.e., N and Fe) varied as for the STD simulation. As expected, the effect on phosphate concentrations (Fig. S9b) and productivity (Fig. S9d) in this sensitivity simulation remains extremely low, with the relative differences to the STD being almost everywhere less than 1%. This overall demonstrates that changes in phosphorus atmospheric deposition do not play a significant role in marine productivity from preindustrial to future periods. 5

Global patterns of productivity change for preindustrial atmospheric deposition inputs
Model calculations demonstrate three major oceanic regions where the reductions in productivity are significantly stronger in PAST compared to PRESENT, i.e., the subtropical gyres of the northern hemisphere Pacific, the Atlantic Ocean, and the northernmost North Pacific (Fig. 3d). The subtropical gyres, however, are the most sensitive to changes in nitrogen input and clearly show the strongest productivity reduction compared to PRESENT. Indeed, the nitrogen concentrations in the subtropics 10 are not much affected (Fig. 2a). This is because good light conditions and warm waters persistently maintain high rates of nutrient consumption, so nitrogen concentrations are already very low in PRESENT. Thus, a change in external nutrient supply feeds immediately into productivity without a significant imprint on nitrogen concentrations.
In the northernmost Pacific, the strong productivity decline in PAST (Fig. 3d) is primarily related to the lowered availability of iron (Fig. 2d), although the reductions in iron deposition remain below 20% (Fig. 1d). However, besides the light conditions, 15 iron availability is the most important factor for limiting productivity in this region (Fig. 4c) compared to nitrogen and phosphorus (Fig. 4a). As a consequence, slight changes in iron supply have a strong impact.

Changes in phytoplankton composition
In the high latitudes, a large part of productivity is related to siliceous diatoms (e.g., Malviya et al., 2016;Uitz et al., 2010), which is accounted for in the model by the low nanophytoplankton to diatoms ratios (Fig. 6b). Accordingly, the overwhelming 20 part of productivity reduction in the northernmost Pacific (Fig. 3d) is related to the decline of diatoms. This is well reflected by the increase of the nanophytoplankton to diatoms ratio for PAST relative to PRESENT (Fig. 6a). In turn, this leads to enhanced silicate concentrations in the North Pacific (Fig. 7a). Part of the unutilized silicate is advected southward via the North Pacific Current and the California Current, leading also to elevated concentrations along the western coast of North America (Fig. 7a). Note, however, that Si atmospheric inputs are here solely depend on the dust deposition fluxes, and thus 25 they have no interannual variability in the model.
A further consequence of the strongly diminished productivity is an accumulation of nitrogen in the subpolar gyre of the North Pacific (Fig. 2a). The nitrogen anomaly is strongest in the southwestern area of the gyre and part of the excess nitrogen is injected into the northern California Current. As a result, a strong positive and wedge-shaped productivity anomaly develops in front of western Canada in PAST (Fig. 3d). This positive anomaly is caused by the increased production of 30 nanophytoplankton productivity (not shown) which dominates in this region, as indicated by higher nanophytoplankton to diatom ratios (Fig. 6b); i.e., north of the wedge lowered iron limits productivity, while south of the wedge nitrate is limiting.
Altogether, this reflects a slight shift from diatom to nanophytoplankton production in the eastern Pacific (north of 40 °N), as indicated by a decline of ~10% of the nanophytoplankton to diatoms concentrations in the upper 20 m (Fig. 7b).
Apart from the northernmost Pacific, the decline in diatom production leads almost everywhere to slightly increased silicate concentrations in PAST (Fig. 7a). Productivity changes in the Southern Ocean, however, remain low for PAST (Fig. 3d). The reason for this is the strong light limitation around Antarctica (Fig. 4b) and the deep mixed layer which suppresses productivity 5 and subsequently builds up a large pool of unutilized nutrients. Part of the unutilized nutrients is advected further north into the Southern Ocean, driving productivity there. Accordingly, the reduced deposition of nitrogen and iron in this area (Figs. 1a,d) have only a slight impact on productivity. Consequently, this region is relatively robust against external nutrient input maintaining stable productivity. A similar effect is seen for the North Atlantic where vigorous exchange with Arctic waters takes place across the Norwegian-Greenland Sea. By contrast, in the subpolar North Pacific, the import of unutilized nutrients 10 from the Arctic is hampered, as the water exchange with polar waters is limited by the shallow Bering Strait and the Aleutian Arc. Therefore, the North Pacific appears the most sensitive to external nutrient inputs compared to other oceanic regions.

Global patterns of productivity change for future atmospheric deposition inputs
For most of the world ocean, productivity changes in FUTURE are in qualitative agreement with PAST, but less pronounced (Figs. 3d,f). This is mostly because both FUTURE and PAST experiments reflect reduced anthropogenic emissions and thus 15 the same mechanisms are involved. The only notable exception compared to PAST (Fig. 3d) is demonstrated in the eastern North Pacific where a strong negative wedge-shaped anomaly is seen in FUTURE (Fig. 3f). This opposite response is related to different iron inputs to the North Pacific (Figs 1d,f). In FUTURE, the reduction of iron atmospheric inputs (Fig. 1f) is by far less strong and in the NE Pacific is even higher than today, thus the productivity increases in the NE Pacific subtropical gyre (Fig. 3f). As a result, more nitrogen is consumed in the subpolar gyre in FUTURE and no nitrogen accumulation takes 20 place as in PAST (Figs. 2a,c). Accordingly, a strong negative nitrogen anomaly develops in the western North Pacific and nitrogen depleted waters are advected southward along the California Current (i.e., opposite to PAST). Altogether, these results imply an extreme sensitivity of the North Pacific against changes in atmospheric iron input. By contrast, the North Atlantic, which is less affected by iron limitation, reflects a widespread decline in productivity mainly controlled by the reduced nitrogen inputs. 25

Biogeochemistry responses to atmospheric organic nutrients
Most marine biogeochemistry studies mainly account for the inorganic fraction as the most important pool of nutrients from the atmospheric pathway. On the other hand, state-of-the-art atmospheric chemistry models nowadays not only efficiently calculate the total dissolved nutrient atmospheric deposition fluxes, but they include the organic part as well, which turns out to be rather important for the total magnitude of the atmospheric input to the ocean (Fig. S3). However, great uncertainty still here the inorganic and organic fractions of N, Fe, and P deposition fluxes to investigate the role of their organic components in marine biogeochemistry. The differences in nitrogen fixation and primary production rates between the STD and ORG simulations are presented in Fig. 8. Note that as for the riverine organic fractions in the model (see Aumont et al., 2015), we also assume here an instant transformation of the atmospheric dissolved organic nitrogen (DON) and organic phosphorus (DOP) inputs to the respective inorganic fractions in the water column. 5 When the organic fraction of the atmospheric nutrients is considered in the model, a modest decrease in the global nitrogen fixation rates of ~0.5 Tg-N yr -1 is calculated for present-day conditions ( Table 1). The increased soluble Fe inputs (~25%) relative to the STD simulations -although smaller compared to relative increases of N (~45%) and P (~50%) -tend to reduce the Fe-limitation in diazotrophs. Consequently, the reduced Fe-limitations make the extra atmospheric inputs of ON to the ocean more effective, overall decreasing the global nitrogen fixation rates in the model (Fig. 8a). Note, however, that the 10 nitrogen fixation is a rather energy-expensive process that is known to be inhibited in the excess of ammonium, in particular.
In the tropical Pacific Ocean, the nitrogen fixation rates for the ORG simulation are significantly more intense compared to STD (up to ~90%) but suppressed (up to ~40%) elsewhere (Fig. 8b). For example, in the Indonesian throughflow and the eastern tropical Atlantic along central Africa, significantly reduced rates (more than 90%) are calculated for the ORG simulation. The lowered nitrogen fixation rates in these regions are mainly due to the additional ON deposition into the ocean. 15 Nitrogen fixation is also decreased in the tropical Atlantic Ocean. On the other hand, increased soluble Fe inputs to the tropical Pacific (Fig. S3d), partially lower the Fe limitation of diazotrophs and increase nitrogen fixation in these remote oceanic regions. Overall, compared to the STD, the present-day global nitrogen fixation rate for the ORG simulation leads to a net decrease by ~0.4% (Table 1).
Primary production increases almost in all ocean basins for the ORG simulation (Fig. 8d), except some parts of the Subpolar 20 Pacific Ocean. In particular, higher rates are calculated in the subpolar Atlantic Ocean (up to 15%). In the N-limited oceanic regions, the increased atmospheric nitrogen deposition (Fig. S3b) directly increases the production rates (Fig. 8d). Such a case is the western subtropical North Pacific, where atmospheric N deposition supports an extra production of up to 15%. The production rates are also increased in the subtropical South Pacific and Atlantic Oceans up to nearly 20%. In total, the primary production increased from ~46.7 Pg-C yr -1 for the STD to ~47.8 Pg-C yr -1 for the ORG (Table 1). Figure 8d points out to 25 regions in the Pacific where production decreased. For the North Pacific, this represents the same mechanism as described above for the differences in primary production rates between PAST and PRESENT (Sect. 3.3). For the ORG simulation, the increased iron input in the Pacific subpolar gyre (Fig. S3d) increases the diatom production leading to higher consumption of nitrate (see Fig. S10b). Subsequent transport of nitrogen diminished waters further south causes a decrease in productivity.
The boundary between decreased and increased bands (Figs. S10b,d) matches, however, the sharp transition from iron 30 limitation to nitrogen limitation (Figs. 4a,c). Increased iron input south of the boundary (i.e., where Fe limits) stimulates the production and diminishes nitrogen. On the other hand, the advective mixing of the N-diminished waters with waters further north decreases the productivity north of the boundary (i.e., where N limits). Overall, the result is the calculated dipole pattern in primary production rates as demonstrated in Fig. 8d.
The main focus of this study is to investigate the effect of nutrient deposition on oceanic primary production. Hence, the presented experiments did not account for the impact of future climate change which could interact or may even mask the effect of changed atmospheric deposition fluxes considered here. Consequently, the here found effects are subject to some uncertainties related to the potential interaction with climate change. For example, climate-induced changes in the global wind 5 system may not only alter atmospheric pathways for nutrients but also impact on oceanic up-and down-welling. Thus, shifts in the seasonal position of trade winds will likewise force shifts in the position of open-ocean and coastal upwelling. These regions are usually nutrient-rich and not particularly sensitive to varying atmospheric nutrient inputs. Displacements of these upwelling positions, as a result of climate change, can increase the sensitivity to external nutrient inputs in regions formerly impacted by upwelling. 10 Several studies have demonstrated that the mid to high latitude areas, such as the North Atlantic and the Arctic, will be more stratified in a future warmer climate (Bindoff et al., 2019;Fu et al., 2016;Gröger et al., 2013;Sein et al., 2018;Steinacher et al., 2010), with negative feedback on vertical mixing and marine primary production due to reduced upward transport of nutrients into the photic zone. Accordingly, primary production in these regions will probably be more sensitive to changed atmospheric deposition rates in the future. Our results, overall, imply only marginal effects in polar regions like the Arctic 15 Ocean. This is certainly robust under the present climate when marine productivity is limited by temperature and sea-ice reducing light conditions in these regions. However, there is a large agreement that climate change will be most severe in the high latitudes, with strong increases in the water temperatures and substantially diminished sea ice cover in the Arctic (Collins et al., 2013). Temperature and sea-ice related light limitation will likely become less important in this region and thus, more nutrients will be recycled in the polar region and less exported equatorward. Consequently, changes in atmospheric transport 20 and deposition of the bioavailable nutrients may play a larger role in a future climate, especially under the high-emission climate scenarios. An example can be seen in the high latitudes of the Southern Ocean around Antarctica where the major amount of surplus DFe is deposited in our PAST and FUTURE experiments (Figs. 1d,f). As expected, the additional DFe availability has nearly no effect on productivity (Figs. 3d,f) as convective mixing and extremely low water temperatures maintain sufficient nutrients and support low productivity under the present-day climate. This may, however, change with 25 altered oceanographic conditions under a future warmer climate. In the northernmost Pacific, known as an HNLC region where iron is the limiting factor, the increased supply of DFe clearly stimulates marine productivity in the PRESENT and FUTURE periods compared to PAST. However, this increase in productivity is likely overestimated here, since our experiments lack climate-induced changes in future stratification which would reduce the nutrient supply from the deep ocean.
The impact of atmospheric organic nutrients on the global oceanic productivity turns out as high (~1 Pg-C yr -1 ; Table 1) as the 30 increase in the present-day primary production since preindustrial times when only the inorganic nutrients' supply is accounted for. However, based on the Krishnamurthy et al. (2009) model estimations, an atmospheric pCO2 declined of about 2.2 ppm due to the modern era's iron and nitrogen inputs to the global ocean compared to preindustrial times can be supported. Hence, the here calculated increase in primary production related to the input of organic nutrients could correspond to an additional decrease in atmospheric pCO2 of ~0.2 ppm, respectively. All changes in nutrient deposition fluxes in the model, are solely driven by changes in the anthropogenic and biomass burning emissions, along with the changes in insoluble to soluble conversions rates due to atmospheric processing. Indeed, the atmospheric deposition fields used in this study did not account for any changes in dust (and bioaerosol) emissions. Instead, they were kept constant to the present-day atmosphere (i.e., the 5 year 2010), although several studies suggest that dust fluxes may be sensitive to climate change and the land-use changes (e.g., Ginoux et al., 2012;Mahowald et al., 2010;Prospero and Lamb, 2003), and thus could be an important driver of the atmospheric nutrient cycles.

Summary and conclusions
This study presents the implementation of state-of-the-art monthly mean atmospheric deposition fields in the global 10 biogeochemistry model PISCES. The model is here run in offline modus, forced by dynamical physical outputs from the physical ocean model NEMO. The newly coupled atmospheric deposition fields considered for this work are all calculated based on a detailed representation of emissions of natural and combustion nutrient-containing aerosols, detailed atmospheric gas-and aqueous-phase chemical schemes, and mineral dissolution processes due to atmospheric acidity and organic ligands.
Another feature tested in the present study is the contribution of organic components to the atmospheric inputs to the global 15 ocean. Moreover, to effectively isolate here the impact of atmospheric deposition on the marine biogeochemistry parameters, the atmospheric CO2 mixing ratio is set to the preindustrial values for all simulations.
For the present day, ~40 Tg-N yr -1 , ~0.28 Tg-Fe yr -1 , and ~0.10 Tg-P yr -1 of inorganic nitrogen, iron, and phosphorus atmospheric inputs to the global ocean are considered in PISCES. This results in a global nitrogen-fixation rate of ~112 Tg-N yr -1 and an integrated primary production of roughly 47 Pg-C yr -1 . Compared to present-day conditions, the lower preindustrial 20 atmospheric nutrient inputs to the ocean result in a weakened primary production of ~2% globally. The decrease in oceanic productivity is supported by the preindustrial decrease in the soluble iron inputs, owing to changes in combustion sources and the atmospheric processing of mineral aerosols, along with the substantial decrease in atmospheric anthropogenic nitrogen inputs. The projected changes in air pollutants under the RCP8.5 emission scenario also result in a modest decrease in marine productivity compared to modern times. Global nitrogen-fixation rates present here a marginal variability, although some 25 notable decreases are calculated for the modern subtropical Pacific and Atlantic gyres.
This work asserts the importance of an explicit representation of the atmospheric nutrients in the context of biogeochemistry modeling, providing also a first assessment of the contribution of another source of atmospheric nutrients than inorganics and, thus, highlights the potential importance of organic nutrients on oceanic productivity. Overall, our main conclusions can be summarized as: 30 1) A general low impact of atmospheric nutrient deposition scenarios on the total marine primary production on a global scale. This is because much of modern productivity is driven by nutrients already recycled in the euphotic zone or by nutrient import from the deep ocean (such as in upwelling regions). Additionally, atmospheric transport appears rather important, as a significant part of nutrient deposition takes place in the northern high latitudes, where light conditions and temperature further limit productivity. Accordingly, even substantial reductions of nitrogen, phosphorus, and iron inputs, ranging between 36 and 51% during the preindustrial period, result in an only modest decline of primary production of about 3% compared to present-day. 5 2) Substantial local productivity changes of up to 20% are found in regions limited by nutrients. The strongest sensitivity against atmospheric nutrients is found for the oligotrophic subtropical gyres of the North Atlantic and the North Pacific, where good light conditions and warm temperatures together with low nutrient concentrations predominate.
Additional atmospheric nutrient input to these regions immediately results in production by increasing the biogenic turnover. 10 3) The North Pacific appears more sensitive to the external nutrient atmospheric deposition compared to other oceanic regions mainly due to two reasons: the strongest deposition changes take place in the northern mid to high latitudes, and that compared to the Southern Ocean and the North Atlantic, the exchange with cold and nutrient-enriched polar waters is limited by land by the shallow Bering Strait and the Aleutian Arc. By contrast, the southern high latitude ocean contains a large amount of unutilized nutrients that are advected further north (to mid-latitudes) making this 15 region more robust against changes in external nutrient input. In agreement, however, with observational evidence from WOA, PISCES exhibits a widespread surplus of nitrogen compared to phosphorus and with respect to the Redfield ratio. Therefore, the applied changes in phosphorus inputs have nearly no impact on primary production in the model. This applies even to the warm water regions, where reductions in atmospheric iron supply limit nitrogen fixation by diazotrophs in both PAST and FUTURE periods. 20 4) The North Pacific turns out to be the most sensitive oceanic to iron atmospheric deposition changes. For the preindustrial period, the lowered input of iron to this region leads to a strong decline of siliceous diatom production leading to an enrichment of silicate, nitrogen, and phosphorus. In turn, this leads to enhanced equatorward transport of nutrients, resulting in elevated production rates of calcareous nanophytoplankton further southeast. 5) Finally, the effect of atmospheric organic nutrient deposition fluxes on the global primary production is calculated 25 roughly as strong as the effect of the present-day increased emissions and atmospheric processing on the oceanic biogeochemistry since preindustrial times when only the inorganic fraction is considered in the model (~1 Pg-C yr -1 ).
Although the overall impact of atmospheric organic nutrient deposition on a global scale is rather low, some stronger changes in regional oceanic productivity are clearly demonstrated in the oligotrophic subtropical gyres. Data availability. Atmospheric nutrient deposition data used for this study are available at Zenodo (doi: 10.5281/zenodo.4017209).
Author contribution. SM prepared the atmospheric input fields, performed the simulations, and the model evaluation. MG prepared the model forcing data. SM and MG wrote the manuscript. All authors contributed to the manuscript preparation.