Simulated Arctic Ocean Response to Doubling of Riverine Carbon and Nutrient Delivery

The Arctic Ocean, more than any other ocean, is influenced by riverine input of carbon and nutrients. That riverine delivery is likely to change with climate change as runoff increases, permafrost thaws, and tree lines advance. But it is unknown to what extent these changes in riverine delivery will affect Arctic Ocean primary production, air‐to‐sea CO2 fluxes, and acidification. To test their sensitivity to changing riverine delivery, we made sensitivity tests using an ocean circulation model coupled to an ocean biogeochemical model. In separate idealized simulations, riverine inputs of dissolved inorganic carbon (CT), dissolved organic carbon (DOC), and nutrients were increased by 1%/year until doubling. Doubling riverine nutrient delivery increased primary production by 11% on average across the Arctic basin and by up to 34–35% locally. Doubling riverine DOC delivery resulted in 90% of that added carbon being lost to the atmosphere, partly because it was imposed that once delivered to the ocean, the riverine DOC is instantaneously remineralized to CT. That additional outgassing, when considered alone, reduced the net ingassing of natural CO2 into the Arctic Ocean by 25% while converting the Siberian shelf seas and the Beaufort Sea from net sinks to net sources of carbon to the atmosphere. The remaining 10% of DOC remained in the Arctic Ocean, but having been converted to CT, it enhanced acidification. Conversely, doubling riverine CT increased the Arctic Ocean's average surface pH by 0.02 because riverine total alkalinity delivery increased at the same rate as riverine CT delivery.


Introduction
It is uncertain how river delivery of carbon and nutrients will change and how these changes will affect the coastal and open ocean (Regnier et al., 2013). The largest of these changes will occur in the Arctic Ocean into which 11% of the global river discharge drains  even though it is the world's smallest ocean, representing only 4% of the global ocean area and 1% of its volume (Jakobsson, 2002).
Total Arctic river discharge has continued to increase since the beginning of the last century, for example, with outflow from the six largest Eurasian rivers draining into the Arctic Ocean growing by 7% from 19397% from to 19997% from (McClelland et al., 2004Peterson et al., 2002). Between 1964 and 2000, river discharge from 16 Eurasian rivers that drain into the Arctic Ocean increased by 11% (McClelland et al., 2006). Conversely, the Canadian river discharge into the Arctic Ocean declined by 10% from 1964 to 2003 (Déry & Wood, 2005), although that was reversed during 1989 to 2007, a period that showed a 15% increase in river discharge (Déry et al., 2009). Further increases in Arctic river discharge are expected given projected future increases in precipitation (Peterson et al., 2002). More precisely, a 16-28% increase of freshwater discharge into the Arctic Ocean during the 21 st century is projected by atmosphere-ocean general circulation models forced under the SRES A1, A2, and B1 scenarios (Lawrence & Slater, 2005;Nohara et al., 2006). Such increases would in turn affect Arctic Ocean circulation and biogeochemistry, for example, leading to increased stratification, decreased vertical mixing, decreased nutrient supply from deeper waters, decreased primary production, and enhanced acidification (Carmack et al., 2015).
Primary production and acidification are also affected by riverine delivery of carbon and nutrients (Semiletov et al., 2016;Tremblay et al., 2015). Out of all the dissolved inorganic carbon (C T ) that is delivered to the global ocean by rivers, 13% to 15% is delivered into the Arctic Ocean . In addition, the ocean delivery of dissolved organic carbon (DOC) from the six largest Arctic rivers is 2.5 times Global Biogeochemical Cycles 10.1029/2019GB006200 larger than that from temperate rivers having similar watershed size and water discharge (Raymond et al., 2007). Riverine carbon and nutrient delivery influences Arctic Ocean biogeochemistry in multiple ways, for example, by increasing primary production due to riverine nutrient delivery (Letscher et al., 2013;Le Fouest et al., 2013, reducing CO 2 uptake over the Siberian shelf seas due to their large riverine DOC delivery Manizza et al., 2011), and enhancing coastal ocean acidification also due to DOC delivery (Semiletov et al., 2016).
Despite the substantial riverine delivery of carbon and nutrients to the Arctic Ocean and the associated potentially large impacts on its biogeochemistry, long-term measurements of ammonium and dissolved organic nitrogen (DON) fluxes in Russian rivers were found to be unreliable (Holmes et al., , 2001. To establish a better database for riverine fluxes to the Arctic Ocean, the Pan-Arctic River Transport of Nutrients, Organic Matter, and Suspended Sediments project (PARTNERS) was launched in 2003. PARTNERS made time-coordinated measurements of nutrients and carbon throughout the year in the six largest Arctic rivers (Ob, Yenisei, Lena, Kolyma, Yukon, and Mackenzie), which cover a combined watershed area of 11.3 × 10 6 km 2 (55% of Arctic watersheds). In 2008, the measurement program of the PARTNERS project was continued via the Arctic Great Rivers Observatory (ArcticGRO).
These observational programs have advanced the ability to assess present-day Arctic riverine fluxes, but of course they do not tell us how those fluxes could change in the future nor the corresponding effects on the biogeochemistry of the Arctic Ocean. Future riverine fluxes of carbon and nutrients are likely to increase due to advancing tree lines (Harsch et al., 2009), deepening of the active permafrost layer (Oelke et al., 2004), and degrading permafrost. Near-surface permafrost is projected to decline from 10.5 to 1.0×10 6 km 2 by 2100 in a fully coupled global climate model (CCSM3) forced under the SRES A2 emission scenario (Lawrence & Slater, 2005), a trend that is expected to enhance riverine delivery of C T (Tank, Frey, et al., 2012;Walvoord & Striegl, 2007) and total alkalinity (A T ) (Drake et al., 2018) to the Arctic Ocean. Simultaneously, there may be an associated 29-46% increase in DOC flux from peatlands, namely, from the West Siberian watersheds, based on the observed relationship between atmospheric temperature and riverine DOC concentrations and projected temperature increases under both the SRES A2 and B2 scenarios (Frey & Smith, 2005). Also projected for the same watersheds are simultaneous increases in concentrations of DON (32-53%), total dissolved nitrogen (30-50%), and total dissolved phosphorus (29-47%; Frey et al., 2007). Nevertheless, Frey and McClelland (2008) question these projected increases in inorganic nitrogen and organic matter delivery because of the large uncertainties associated with river discharge projections. Faced with uncertainties, scientists have used idealized forcing scenarios to characterize the response of complex systems. For example, Friedlingstein et al. (2001) used idealized simulations to describe the response of the land and ocean sinks to increasing atmospheric CO 2 and changing climate. That is, with a system that combines a coupled ocean-atmosphere general circulation model and models of the carbon cycle on land and in the ocean, they increased atmospheric CO 2 by 1%/year and assessed the amount of atmospheric CO 2 that was taken up by land and by ocean. The responses of these land and ocean sinks were then divided into those from increasing atmospheric CO 2 and from changing climate by making two different simulations. Likewise, Boer and Arora (2010) as well as Roy et al. (2011) calculated the same responses but with more plausible socioeconomic scenarios. Yet calculated climate sensitivities depend on the scenario, so there has been a return to using the classical 1%/year idealized increase scenario (Arora et al., 2013).
This same idealized approach could be transposed to assess sensitivities of how changes in atmospheric CO 2 and riverine input of carbon and nutrients, driven in part by climate change, will alter projected changes in ocean acidification (Steinacher et al., 2009;Steiner et al., 2013), air-to-sea CO 2 fluxes (Bates et al., 2006), ocean C T (Anderson & Kaltin, 2001), and primary production (Vancoppenolle et al., 2013), all of which may directly or indirectly affect the marine ecosystem (Darnis et al., 2012;Riebesell et al., 2013). Quantifying such sensitivities would offer common ground from which to compare models and gauge their developments regarding how river fluxes affect Arctic Ocean biogeochemistry.
For example, models disagree on the sign of the change in primary production during the 21st century (Vancoppenolle et al., 2013). They also disagree among themselves concerning the extent of future shoaling of the deep aragonite saturation horizon . Likewise, the sign of the future change in the air-to-sea CO 2 flux in the Arctic Ocean projected by models (Roy et al., 2011) is opposite that estimated by one observational study (Cai et al., 2010), while it agrees with another (Bates et al., 2006). Yet none of the model studies consider increases in the partial pressure of CO 2 (pCO 2 ) from changes in riverine DOC. Few earth Global Biogeochemical Cycles 10.1029/2019GB006200 Figure 1. Arctic regional seas and its six major rivers, out of which all but the Yukon drain into the Arctic Ocean. system models account for riverine input of carbon and nutrients. Those that do only do so in a rudimentary way even though observations indicate that such riverine input affects Arctic Ocean acidification (Chierici & Fransson, 2009;Semiletov et al., 2016), air-to-sea CO 2 fluxes (Cai et al., 2010;Manizza et al., 2011), and primary production (Le Fouest et al., 2015;Tank, Manizza, et al., 2012).
Our aim here is to assess how changes in riverine fluxes affect simulated biogeochemistry, in part by providing sensitivities that quantify the effects from changes in riverine delivery of carbon and nutrients and those from increasing atmospheric CO 2 . For simplicity, this study neglects effects of climate change on Arctic Ocean biogeochemistry, such as warming, stratification, and sea ice melt.

Arctic Regions
To assess spatial patterns in how changes in riverine carbon and nutrient delivery affect its biogeochemistry, the Arctic Ocean is divided into eight regions (Figure 1). The Central Arctic includes all waters where the seafloor is deeper than 500 m. The remaining area is then divided into seven coastal seas, classed either as "exterior" (Barents Sea, Chukchi Sea, and Canadian Arctic Archipelago [CAA]) because of their direct exchange with the Atlantic or Pacific Ocean or "interior" (Kara Sea, Laptev Sea, East Siberian Sea, and Beaufort Sea) because they have no such exchange.

Coupled Ocean-Biogeochemical Model
This study relies on the ocean general circulation modeling platform Nucleus for European Modelling of the Ocean (NEMO), more specifically its version "v3.6 stable." Here we use the full NEMO system, consisting of the ocean general circulation model OPA (Madec, 2008), the Louvain-la-Neuve Sea Ice Model (LIM3.6) (Rousset et al., 2015), and the "Tracers in the Ocean Paradigm" (TOP) model. In our case, TOP couples NEMO-LIM to the biogeochemical model "Pelagic Interactions Scheme for Carbon and Ecosystem Studies" (PISCES-v2; Aumont et al., 2015).
Simulations were made with the same global configuration of NEMO known as ORCA1 (1 • nominal horizontal resolution), having a normal Mercator grid south of 20 • N but a distorted grid north of that boundary Global Biogeochemical Cycles 10.1029/2019GB006200 to avoid the standard grid singularity at the North Pole (over ocean). Instead, that North Pole singularity is replaced for numerical efficiency by two grid singularities over land (over North America and over Eurasia; Madec & Imbard, 1996). That distortion also causes the model's horizontal resolution to be higher in the Arctic Ocean, where the horizontal grid length varies from 25 to 63 km. Vertically, the model is split into 75-depth levels whose thicknesses increase with depth from 1 m (level 1) to 204 m (level 74). The depth of the deepest cell (level 74) can reach up to 408 m, being extended into level 75 as a function of the bathymetry (partial steps; Barnier et al., 2006). Also, at other vertical levels, the partial steps approach allows the depth of the deepest grid cell to be variable and thus permits a better representation of the ocean bathymetry. The ORCA1 global bathymetry map is derived from three different sources: (1) the 2-min ETOPO2 bathymetry map from the National Geophysical Data Center, applied over most of the ocean (Smith & Sandwell, 1997), (2) the IBCAO bathymetric data, applied in the Arctic (Jakobsson et al., 2000), and (3) the BEDMAP bathymetric data, applied south of 72 • S (Lythe & Vaughan, 2001). Because ORCA1 does not explicitly resolve ocean eddies, subgrid-scale eddy effects are parameterized by implementing the Gent and McWilliams (1990) scheme with an eddy diffusion coefficient of 1,000 m 2 /s. The corresponding lateral diffusivity coefficient is 1,000 m 2 /s, while its lateral viscosity coefficient is 2 × 10 4 m 2 /s. A detailed description of the biogeochemical model PISCES is provided by Aumont et al. (2015). Briefly, it simulates four plankton types (nanophytoplankton, diatoms, microzooplankton, and mesozooplankton) as well as the biogeochemical cycles of the main nutrients (N, P, Fe, and Si), C T , A T , and dissolved O 2 . Its total net primary production (NPP) depends on temperature and is limited by light and nutrients. In PISCES, N, P, and Fe limit growth of all phytoplankton, while Si also limits growth of diatoms. For all plankton, the C:N:P molar ratio of organic matter is held constant at 122:16:1 (Takahashi et al., 1985), while the O 2 :C molar ratio is fixed at 1.34 (Körtzinger et al., 2001). The same C:N:P ratio is also fixed for the PISCES nonliving compartments of marine semilabile dissolved organic matter as well as small and large sinking particles. The PISCES model explicitly simulates, as separate tracers, concentrations of phytoplankton chlorophyll, Fe, and Si (for diatoms only). The air-to-sea CO 2 flux in PISCES is calculated from the air-sea difference in the partial pressure of CO 2 , wind speed, sea ice fraction, and CO 2 solubility and Schmidt number as summarized by Bourgeois et al. (2016). In PISCES, calcite is the only form of CaCO 3 that is explicitly simulated. It is transported as a passive tracer in the model and its internal sources and sinks are dissolution and precipitation. That "CaCO 3 concentration" is not used to compute the aragonite saturation state (Ω arag ), which would be erroneous. Rather, Ω arag is calculated offline from [Ca 2+ ] and [CO 2− 3 ] concentrations using the routines from Orr and Epitalon (2015) with equilibrium constants recommended for best practices and simulated temperature, salinity, C T , A T , total dissolved silicon (Si T ), and total dissolved inorganic phosphorus (P T ).

River Input
At the beginning of all simulations, the river inputs into the Arctic Ocean are based on annual fluxes of terrigenous DOC, dissolved inorganic nitrogen (DIN), DON, P T , dissolved organic phosphate (DOP), and Si T from the Global NEWS 2 model (GN2; Mayorga et al., 2010) and C T from the Global Erosion Model (GEM; Ludwig et al., 1998). The GN2 model is a composite of independent submodels for dissolved inorganic, dissolved organic, and particulate C, N, and P, as well as dissolved Si. The GN2 submodels for dissolved elements use a unified formulation. All GN2 submodels use hydrological and physical factors, hydrography, and basin characteristics to estimate annual riverine exports of DOC and nutrients. The GN2 submodels consider both natural processes and anthropogenic activities. GEM establishes a relationship between C T yield and hydroclimatic and geomorphological factors for four different climatic zones in 60 different river basins. These relationships depend on the climatic zone and are then applied to all rivers around the world. In our simulations, lateral boundary conditions at river mouths are applied based on contemporary annual river fluxes from GEM's C T (R C T ) and GN2's DOC (R DOC ), DIN and DON (R N ), P T and DOP (R P ), and Si (R Si ). Along with the corresponding freshwater discharge rates given with GEM and GN2, riverine concentrations of carbon and nutrients are calculated. These concentrations are then multiplied with the monthly river discharge from Dai and Trenberth (2002) used in NEMO. The resulting river fluxes are presented in Table 1.
For simplicity, in the ORCA1-PISCES model it is assumed that the riverine A T :C T ratio is 1.0. Conversely, the observed A T :C T ratio in rivers ranges from 0.6 in the Congo River (Wang et al., 2013) to 1.1 in the Delaware Estuary (Joesoef et al., 2017). In the Mississippi River, that ratio is 1.0 (Cai, 2003). For the six largest Arctic rivers,  observed an average carbonate alkalinity (A C ):C T ratio of 0.91 (0.72-0.94). Although they measured A T , they report A C , after subtracting the measured contribution from organic acids, assuming that alkalinity contributions from P T and Si T are negligible. Thus, their Table 1 River-to-Ocean Fluxes of C T Derived From  Estimates of R C T , R DOC , R N , R P , and R Si Manizza et al., 2009;  Note. Fluxes into the Arctic Ocean and its regional seas are given as the cumulative amount draining into each region from all Arctic rivers. GEM = Global Erosion Model; DOC = dissolved organic carbon; CAA = Canadian Arctic Archipelago.
A C is directly comparable to the simulated A T in PISCES, which neglects organic acids. By using a globally constant A T :C T ratio of 1.0, we overestimate the riverine A T flux (R A T ) by 6-28% depending on the river.
PISCES also models the influence of R DOC simplistically. In the standard version of PISCES, the imposed river flux of terrigenous DOC (R DOC ) is assumed to be extremely labile, being remineralized instantaneously, that is, it is added to the ocean in the form of already remineralized carbon (C T ; Aumont et al., 2015). That simplicity is maintained here, considering that as a first sensitivity assessment, our preference is to provide a limit rather than a best estimate for ocean behavior. Indeed, there are large uncertainties concerning the lability of terrigenous DOC. In North American rivers, Holmes et al. (2008) estimate that 20-40% of terrigenous DOC remineralizes within 3 months, whereas in Eurasian rivers Kaiser et al. (2017) estimate that close to 50% of terrigenous DOC remineralizes within a year. That partitioning may also change in the near future, as more old, labile terrestrial carbon becomes available due to thawing of permafrost (Vonk et al., 2013). Similarly, riverine fluxes of DON and DOP were added to the ocean as DIN and P T .
Given our two simplifications, that is, that the riverine flux of A T (R A T ) equals R C T and that R DOC is instantaneously converted to C T as it is added to the ocean, the effective ratio R ⋆ , imposed as the river flux boundary condition for PISCES, is where are the river flux boundary conditions for PISCES calculated from the river fluxes derived from GEM (R C T ) and Global NEWS 2 (R DOC ). Thus, an increase in while an increase in R DOC would reduce it. To illustrate the uncertainties associated with the two simplifications, we compared our calculated R ⋆ boundary condition for the five major rivers that drain directly into the Arctic Ocean to what it would have been had we used a more realistic R DOC lability of 50% (Kaiser et al., 2017) and the R A T :R C T from , which varies among those rivers ( Figure S1 in the supporting information). Relative to our simplified case, an R DOC lability of 50% would increase the R ⋆

A T
:R ⋆ C T ratio, while imposing a R A T :R C T below one, as observed, would reduce it. Combined, the two simplifications cause the R ⋆ A T :R ⋆ C T ratio to be overestimated in the CTL simulation by 6% in the Kolyma River and to be underestimated by 5-22% in the remaining four rivers. Due to our two simplifications, the change ratio when riverine DOC is doubled is underestimated by 7% for the Ob and overestimated by 26-55% in the remaining four rivers. However, our two simplifications do not alter the direction of the change in the R ⋆ ratio, as discussed in section 4.2. Invariably, that ratio declines when the riverine DOC flux is doubled, and it increases when the riverine C T flux is doubled.

Transient Simulations
We made five 75-year simulations, each with the same physical forcing, the daily climatological DRAKKAR Forcing Set 4.4 (DFS4.4) (Brodeau et al., 2010). The DFS4.4 forcing includes historical reanalyses of atmospheric air temperature and humidity at 2 m, zonal and meridional wind fields at 10 m, downward shortwave and longwave radiation at 2 m, and the net surface freshwater flux (precipitation minus evaporation). The first version of DFS4.4 is based on the ERA40 reanalysis (Uppala et al., 2005) and covers 45 years . DFS4.4 was extended until 2012 using ERA-interim reanalysis (Dee et al., 2011) thus covering 55 years in total. For our 75-year simulations, after the first 55 years we reused the initial 20 years of DFS4.4. This relooping of DFS4.4 is facilitated by efforts made during its construction to make adjustments for global and regional biases and to reduce time discontinuities in data from different sources, which would otherwise result in spurious trends (Brodeau et al., 2010).
To quantify the effect of changes in riverine delivery on Arctic biogeochemistry, we separately altered one aspect of each of the simulations: (1) a preindustrial control simulation (CTL) with a constant atmospheric CO 2 fixed at 284.7 ppm (1850) having constant riverine input as prescribed in section 3.1; (2) a simulation just like CTL except that atmospheric CO 2 is increased by 1%/year (CO 2 ); (3) a simulation just like CO 2 except that R DOC is also increased by 1%/year (ROC), which is instantaneously and completely remineralized as it enters the ocean, thus increasing R ⋆ C T ; (4) a simulation just like ROC except that R C T is also increased by 1%/year (RIC), which increases R ⋆ A T as well as R ⋆ C T ; and (5) a simulation just like RIC except that the nutrient river flux (R NUT ) is also increased by 1%/year (RUT; Table 2). The 1%/year increase in riverine inputs concerns all rivers, not just those in the Arctic. That rate of increase leads to a doubling after 70 years. Taking differences between simulations allows us to distinguish the effects of increases in atmospheric CO 2 , R DOC , R C T , and R NUT on Arctic Ocean biogeochemistry. These differences between simulations are referred to as ΔCO 2 (CO 2 -CTL), ΔROC (ROC-CO 2 ), ΔRIC (RIC-ROC), and ΔRUT (RUT-RIC).
To smooth out interannual variations, all model results were compared in terms of their averages over simulated years 66-75, that is, centered around the time when the simulated atmospheric CO 2 doubled from the imposed 1%/year increase. To assess air-to-sea CO 2 fluxes, NPP, pH, and Ω arag , we also calculated 10-year averages for years 26 to 35 when atmospheric CO 2 varies between 369 and 399 ppm, similar to values observed at the beginning of this century. Initial conditions for atmospheric CO 2 represent the preindustrial level, while those for river input and climate are based on modern data.
Simulated biogeochemical processes, such as the air-to-sea CO 2 flux and NPP, are affected by sea ice. Previous evaluation of the sea ice model component LIM3 coupled to the same ORCA1 configuration and forced by an updated version of the DRAKKAR Forcing Set 5.2 indicates that the simulated mean seasonal cycle and interannual variations of sea ice extent in the Arctic Ocean are much like those observed (Uotila et al., 2017). This simulated maximum sea ice extent in winter is ∼7% larger than observed while minimum sea ice extent in summer generally agrees with observations. Carbon storage Air-to-sea CO 2 flux NPP pH Note. Numbers in parentheses indicate the equations where these factors are defined in the text. NPP = net primary production; DOC = dissolved organic carbon.

Sensitivity Factors
Following the approach first proposed by Friedlingstein et al. (2003), we assessed the sensitivities of land and ocean carbon sinks to increasing atmospheric CO 2 and to anthropogenic climate change. This approach was extended here to calculate sensitivities of the Arctic Ocean carbon sink to the increasing riverine delivery of carbon and nutrients from rivers draining directly into the Arctic Ocean (R DOC , R C T , and R NUT ; Table 3). Furthermore, these sensitivities were not only calculated for the ocean carbon storage and air-to-sea CO 2 flux but also for net primary production and pH.
To quantify the sensitivity of global ocean carbon storage to the atmospheric CO 2 increase ( global ΔCO 2 ), our approach exactly follows that of Friedlingstein et al. (2006) global where ΔC A is the change in the atmospheric CO 2 mixing ratio (ppm) between the beginning and end of the simulation, ΔC global O is the change of C T in the global ocean during the same period, and ΔF global is the difference in global air-to-sea CO 2 flux in response to increasing atmospheric CO 2 . At the global scale, the time-integrated air-to-sea CO 2 flux and the change in carbon storage are identical. Regionally though, they are unlikely to be the same. In the Arctic Ocean, much of the anthropogenic carbon enters laterally (Terhaar et al., 2019a), causing these two diagnostics to differ. To distinguish between them in the Arctic Ocean, is used to indicate the sensitivity of the air-to-sea CO 2 flux (equation (3)), while is used for the sensitivity of ocean carbon storage where ΔF CO 2 is the difference in air-to-sea CO 2 flux in response to increasing atmospheric CO 2 areally integrated over the Arctic Ocean and ΔC is the change in C T storage within the Arctic Ocean in response to an increase in atmospheric CO 2 .
Likewise, the sensitivities of the air-to-sea CO 2 flux to changes in R DOC , R C T , and R NUT are as follows: where ΔR DOC , ΔR C T , and ΔR NUT are differences of riverine input of each of those species into the Arctic Ocean between a given year and that at the beginning of each simulation, while ΔF ROC , ΔF RIC , and ΔF RUT are the differences in air-to-sea CO 2 fluxes in response to increases in R DOC , R C T , and R NUT between the same two times after integrating areally over the Arctic Ocean. Similarly, the sensitivities of Arctic Ocean carbon storage to increasing riverine inputs of DOC, C T , and nutrients are defined as are the time differences in Arctic Ocean C T , integrated vertically and areally over the Arctic, induced by the corresponding increases in riverine input. The ΔRUT does not account for ocean POC and organic carbon burial.
It follows that the sensitivity of NPP ( ) to changes in riverine nutrients is where ΔNPP RUT is the change in the vertically and horizontally integrated primary production from increased riverine nutrients. In the same way, sensitivities of pH ( ) to changes in atmospheric CO 2 and inputs of riverine DOC and C T are where ΔpH CO 2 , ΔpH ROC , and ΔpH RIC are the areally averaged changes in surface pH in response to increases in atmospheric CO 2 , river input of DOC, and river input of C T , respectively. To calculate these basinwide differences, the pH values at every grid cell in the Arctic Ocean were first converted to hydrogen ion concentrations [H + ] on the total scale. The [H + ] was then areally averaged at each time; the results were converted back to pH units, and the difference between the two simulations was taken.

Riverine Forcing
To evaluate the Arctic river flux estimates derived from GN2 and GEM, they were compared to data-based fluxes from the five largest Arctic rivers that drain directly into the Arctic Ocean (Table 1), that is, R C T is from , R DOC from Holmes et al. (2012) and Manizza et al. (2009), and R NUT from Holmes et al. (2012). The model's GEM-derived total R C T flux to our Arctic Ocean domain (50 Tg C/year) falls within the uncertainty range of 's data-based estimate (41 ± 10 Tg C/year) from their extrapolation of the fluxes from the six major Arctic rivers to all rivers that discharge into the Arctic Ocean. The GEM-derived R C T fluxes typically overestimate the extrapolated fluxes but by no more than 12% in the East Siberian Sea, East Siberian Sea, the Chukchi Sea, the Beaufort Sea, and the CAA. Larger overestimates are found for the Barents Sea (+51%) and the Laptev Sea (+126%), while the general tendency is reversed in the Kara Sea (−32%).
The model's GN2-derived total R DOC flux to our Arctic Ocean domain (20.3 Tg C/year) is 32% smaller than the data-based estimate (Manizza et al., 2009). The largest differences are located in the Barents Sea and . Simulated (left column) and data-based (right column) net primary production (NPP; integrated vertically), air-to-sea CO 2 flux, and surface pH and Ω arag (top to bottom). Data-based depth-integrated NPP is derived from remotely sensed ocean color from SeaWiFS over 1998(Arrigo & van Dijken, 2011. The data-based, annual mean air-to-sea fluxes of total CO 2 from 1997 to 2013 are estimates derived with self-organizing pCO 2 maps (Yasunaka et al., 2016). Data-based surface pH and Ω arag are from the GLODAPv2 gridded data product that normalizes observations to 2002 (Lauvset et al., 2016). The modeled air-to-sea CO 2 flux is the decadal average over years 26-35, over which the average atmospheric CO 2 was similar to that during the time span of the data-based estimates. The modeled pH and Ω arag are averaged over July, August, and September to compare with the summer-biased GLODAPv2 data in the Arctic; conversely, the simulated air-to-sea CO 2 fluxes and NPP are annual average conditions as are the data-based estimates. Laptev Sea, where the model's R DOC fluxes are 52-56% smaller than the data-based estimates. The model's underestimated R DOC flux into the Laptev Sea is largely caused by the underestimated GN2-based R DOC flux from the Lena River that is 55% smaller than the data-based estimate. Furthermore, R DOC fluxes for the Ob and Yenisei are 20-37% smaller, while those for the Kolyma agree with the data-based estimates. Unlike for those Siberian rivers, the GN2-based R DOC flux for Canada's Mackenzie River is 36% larger than the data-based estimate.
For R NUT , data-based fluxes have not been extrapolated to all Arctic rivers. They are available only for the five largest rivers draining directly into the Arctic Ocean . For R N , the GN2-based flux for the five observed rivers is 41-169% larger than data-based estimates. For R P , the GN2-based flux estimates for the Yenisei, Lena, Kolyma, and Mackenzie are 50-200% larger than data-based estimates, while that for the Ob is 12% smaller. For R Si , the GN2-based estimates for the Lena and Kolyma Rivers are 15-33% smaller than data-based estimates, while those for the Ob, Yenisei, and Mackenzie Rivers are 18-107% larger.

NPP
Simulated NPP from the CTL simulation was compared to the data product from Hill et al. (2013) derived from remotely sensed ocean color and observed vertical profiles of in situ chlorophyll a. From this data product, the estimated annual mean NPP, integrated vertically and areally over the Arctic Ocean, is 433 Tg C/year. The simulated NPP is much less (165 Tg C/year) but captures similar regional patterns ( Figure 2 and Table 4).
Regionally, the largest model-data differences are found in the three exterior shelf seas: out of the basinwide integrated bias, 49% (132 Tg C/year) is located in the Barents Sea, 31% (84 Tg C/year) in the CAA, and 24% (63 Tg C/year) in the Chukchi Sea (Table 4). The Barents Sea is strongly influenced by inflow from the adjacent Atlantic Ocean while the Chukchi Sea is influenced by inflow from the adjacent Pacific Ocean. Inflow of nutrients from those neighboring oceans fuels at least 20% of the Arctic Ocean NPP (Popova et al., 2013). That nutrient inflow may well be underestimated in the coarse-resolution model used here (ORCA1). Lateral water fluxes into the Arctic Ocean are underestimated using a coarser-resolution version (ORCA2) of the same NEMO-PISCES model, and they are improved when using higher-resolution eddying versions ORCA05 and ORCA025 (Terhaar et al., 2019a). Thus, our weak modeled lateral nutrient inflow may result in the low simulated NPP in the Barents Sea and Chukchi Sea. In the CAA, the discrepancy appears to be mainly caused by differences in the definition of the regional borders of the CAA, which is defined here to be bounded on the south by the Baffin Bay, but for Hill et al. (2013) it extends further south to the Davis Strait. The latter domain includes additional ice-free areas with strong primary production. Unlike for these exterior seas, the integrated simulated NPP in the interior seas is 4% (2 Tg C/year) smaller than the data-based estimates. The total underestimation in these exterior and interior seas adds up to more than 100% of the total bias, because the latter is reduced by an overestimation of simulated NPP in the Central Arctic by 540% (12 Tg C/year). Although most of the difference between simulated and data-based NPP may well be caused by model shortcomings, estimating NPP from remotely sensed ocean color remains challenging in the Arctic Ocean because of the large amount of colored dissolved organic matter, which is difficult to distinguish from chlorophyll a (Matsuoka et al., 2011). Therefore, part of the discrepancy between data-based and simulated NPP could also originate from observational artifacts.

Total Air-to-Sea CO 2 Flux
At present, all subregions of the Arctic Ocean take up both natural and anthropogenic carbon from the atmosphere. Eventually, most of that absorbed carbon is transported out of Arctic Ocean laterally (Bates & Mathis, 2009;Yasunaka et al., 2016), although some local outgassing is observed in the Pacific dominated sectors immediately after sea ice retreat in spring (Arrigo et al., 2010) and on the Siberian shelf in summer after high river discharge of terrigenous DOC ). The annually average total air-to-sea CO 2 flux in the Arctic Ocean is estimated to be 66-199 Tg C/year over (Bates & Mathis, 2009). The simulated Arctic Ocean air-to-sea CO 2 flux at the same level of atmospheric CO 2 is 92 Tg C/year falling within the data-based range (Table 5). Three fourths of that air-to-sea CO 2 flux occurs in the Barents Sea. The next largest simulated air-to-sea CO 2 flux occurs in the Chukchi Sea but is 10 times smaller. The same regional pattern is displayed by data-based estimates with the Barents Sea absorbing 44-77 Tg C/year and the Chukchi Sea taking up 11-53 Tg C/year. This underestimation of the total air-to-sea CO 2 flux in the Chukchi Sea is consistent with the same region's underestimated NPP (Table 4).

Carbonate Chemistry
We evaluate the simulated surface ocean pH by comparing preindustrial, present-day, and future pH estimates to the simulated pH at the corresponding atmospheric CO 2 levels in the CO 2 simulation. Our simulated preindustrial surface pH averaged over the Arctic Ocean is 8.17, which is 0.06 less than that from a previous study with a coupled carbon-climate model (Steinacher et al., 2009).
When doubling atmospheric CO 2 in our model (569 ppm), average surface ocean pH decreases by 0.3. One third of this decrease occurs by the time that the model's atmospheric CO 2 reaches 347-379 ppm (equivalent to observed values over 1986-2005). The 0.1 decrease in pH when atmospheric CO 2 reaches 347-379 ppm agrees with the data-based estimate for that historical surface pH change for the Arctic Ocean (Anderson et al., 2010). The further 0.2 pH reduction by the time of atmospheric CO 2 doubling agrees with the projected decrease for Arctic Ocean surface pH calculated by an ensemble of CMIP5 models forced under the RCP4.5 scenario , where atmospheric CO 2 reached 583 ppm 2100.
The simulated mean Arctic surface pH averaged over years 26-35, when the atmospheric CO 2 forcing averages 374 ppm, was further compared to the gridded GLODAPv2 climatology (Lauvset et al., 2016). Although the simulated surface pH is on average 0.01 too high, the model and data product both indicate notably lower pH in regions with high freshwater input (on the East Siberian shelf, along the transpolar drift, and near the mouths of the Ob and Yenisei Rivers) and higher pH along the east coast of Greenland, in the Kara Sea, and near the St. Anna trough (Figure 2). Besides model deficiencies, differences between the model and the data-based product might be caused by the relatively sparse observations in space and time.
Regional averages of simulated preindustrial Ω arag are all supersaturated (Table 6), consistent with other studies (Anderson et al., 2010;Steinacher et al., 2009), but are spatially heterogeneous, varying from 1.17 to 2.01. The regions corresponding to the largest values of Ω arag are those that are strongly influenced by outside inflow, that is, the Barents Sea (Ω arag = 2.01) adjacent to the Atlantic and the Chukchi Sea (Ω arag = 1.66) adjacent to the Pacific. The regions corresponding to low values of Ω arag are more influenced by riverine input, for example, the Laptev Sea (Ω arag = 1.17) and the East Siberian Sea (Ω arag = 1.17) ( Table 6). The simulated spatial patterns in Ω arag are similar to those found in the GLODAPv2 gridded data product (Lauvset et al., 2016; Figure 2) as well as in several other observational studies, all of which indicate near-zero Ω arag values near the mouths of the Mackenzie (Chierici & Fransson, 2009), Yukon (Mathis et al., 2011), and Lena Rivers (Semiletov et al., 2016).
The low Ω arag in the coastal Arctic Ocean is driven by the dilution from Arctic river waters, with not only relatively low A T :C T ratios (0.7-0.9) but also low A T (500-1,600 μmol/kg) and C T (650-1,700 μmol/kg; . Typical open ocean values are higher for the A T :C T ratio (1.1-1.2) and for both A T (2,300 μmol/kg) and C T (2,150 μmol/kg). Out of the five major rivers that drain into the Arctic Ocean, it is the Kolyma River that exhibits the lowest A T and C T (500-1,000 μmol/kg) concentrations and the lowest A T :C T ratio (0.73), perhaps because of rapid degradation of large amounts of labile terrigenous DOC originating from the Arctic's largest watershed that is entirely covered by permafrost (Mann et al., 2015).  Table 4). The large increase in the Barents Sea represents 27% of total Arctic NPP increase, although the increase in the R NUT boundary condition represents only 10% of the total Arctic R NUT increase (Table 1). This discrepancy suggests that region's increase in NPP is mainly driven by the additional influx of nutrients from outside the Arctic, given that the imposed R NUT boundary condition increases in all rivers across the globe. The large increase in the Kara Sea represents 24% of the total Arctic NPP increase. As opposed to the Barents Sea, the increase of R NUT into the Kara Sea is 47% of the total Arctic R NUT increase. This large amount of extra nutrients thus solely explains Kara Sea's strong increase in NPP. Although the largest absolute changes in NPP occur in the Barents and Kara Seas, the largest relative changes occur in the Laptev and Beaufort Seas (+34% and +35%, respectively; Table 4).

Changes Due to Increasing Riverine Input
Basinwide, the largest absolute changes in NPP are typically found within 100-200 km of the coastline, where they are often twice as large as average changes in adjacent waters further offshore (Figure 3). Yet the coastal-to-open ocean gradients in NPP may well be overestimated because in our simulations it is imposed that DON and DOP river fluxes (R N and R P ) are instantaneously transformed into inorganic N and P as they enter the Arctic Ocean (section 3.1). In the real ocean, portions of the DON and DOP are transported offshore before being transformed into inorganic N and P, thus reducing the coastal-to-open ocean gradients in inorganic nutrients and NPP. Changes in NPP from the doubling of atmospheric CO 2 , R C T , and R DOC are generally more than a hundred times smaller than changes induced by doubling R NUT . Indeed, the version of PISCES used here does not account for any direct effect of increasing C T on phytoplankton productivity. The only indirect effects of changes in C T on phytoplankton productivity are those involving carbonate dissolution and its implications on iron scavenging.
The simulated increase in NPP from the increase in R NUT may be put into context by comparing it to NPP changes from climate-related reductions in sea ice projected by an ensemble of 11 CMIP5 models (Vancoppenolle et al., 2013). When these models were forced under the RCP8.5 scenario, the projected changes in Arctic Ocean NPP ranged from −110 to 253 Tg N/year for the average over 2080-2099 minus that over 1980-1999. Over that period, their average increase (58 Tg C/year) from climate change is around 3 times larger than the Arctic's basinwide increase due to the doubling of nutrient delivery in our RUT simulation. The NPP increase from increased river nutrient fluxes also reduces surface C T and thus surface ocean pCO 2 , consequently enhancing the Arctic Ocean's air-to-sea CO 2 flux. The regions with the largest simulated increase in NPP (Barents Sea, Kara Sea, Laptev Sea, and the CAA; Table 4) also have the largest increase in the air-to-sea CO 2 flux (>0.6 Tg C/year for each) at year 70 (Table 5). These ties between NPP and air-to-sea CO 2 flux are similar to those seen in the other regional seas (Figures 3 and 4). At the time of R NUT doubling, the change in NPP (18 Tg C/year) enhances the air-to-sea CO 2 flux and hence the Arctic Ocean's carbon storage by 1.83 Tg C/year of which 74% is stored in the Central Arctic Ocean, 11% in the CAA, and 6% in the Barents Sea. The remaining 7% is spread over the other Arctic shelf seas.

Effect of Riverine Carbon Increase
The simulated doubling of R DOC reduces the net Arctic Ocean air-to-sea CO 2 flux by 18.2 Tg C/year, compensating 94% of the increase in the Arctic Ocean's air-to-sea CO 2 flux from the doubling of atmospheric CO 2 (Table 5). However, even if the air-to-sea CO 2 flux from the doubling of atmospheric CO 2 was completely compensated by R DOC , the Arctic Ocean would still remain a sink of anthropogenic carbon, which mainly enters the Arctic Ocean laterally from the adjacent ocean basins (Olsen et al., 2015;Terhaar et al., 2019a).
Out of the added R DOC , which is instantaneously converted to C T , 90% is lost from the Arctic Ocean via outgassing of CO 2 , 7% is transported out of the Arctic Ocean laterally, and 3% remains there. That outgassing is strongest in the Kara Sea (6.1 Tg C/year), which changes from a net sink to a net source of CO 2 ( Table 5). The same change from CO 2 sink to source occurs in the East Siberian Sea. Although the Laptev and Beaufort Seas already exhibit outgassing at the end of the CO 2 simulation, the doubling of R DOC increases that outgassing by up to several times. Thus, when R DOC , is doubled, all interior Arctic shelf seas become sources of Global Biogeochemical Cycles 10.1029/2019GB006200 carbon to the atmosphere even when atmospheric CO 2 levels have also doubled. Nevertheless, the Barents Sea, the Chukchi Sea, and the Arctic Ocean as a whole remain as sinks of atmospheric CO 2 .
When R DOC is doubled, the shelf waters immediately adjacent to large river mouths become small hot spots of CO 2 outgassing with intensities that are roughly 2 to 4 times as strong as associated outgassing from each surrounding regional sea (Figure 4). Away from river mouths, the associated outgassing declines rapidly as river waters are diluted in the sea and surface ocean pCO 2 plummets.
The assumption that the R DOC boundary flux for PISCES is instantaneously converted to a flux of C T as it enters the Arctic Ocean is merely a first limiting case for this idealized case study; it is much simpler than what occurs in the real world. From observed concentration gradients of terrigenous DOC in the Beaufort Gyre, Hansell et al. (2004) estimate that the half-life of terrigenous DOC in the Arctic Ocean is 7.1 ± 3.0 year, assuming exponential degradation and average residence time of river waters in the Arctic Ocean of 11 to 15 years based on isotopic water mass tracers. Hansell et al. (2004) further estimate that 21-32% of R DOC is exported laterally to the North Atlantic before being remineralized to C T in the Arctic. More recent studies suggest larger uncertainties, that is, that 20-50% of terrigenous DOC is remineralized in estuaries or the shelf seas Kaiser et al., 2017;Letscher et al., 2011;Spencer et al., 2009). Besides those large uncertainties, the future rate of degradation of terrestrial carbon may increase as old, terrestrial DOC, which is more labile, is mobilized during continuous thawing of permafrost Letscher et al., 2011;Vonk et al., 2013).
The model assumption that R DOC is added to PISCES as C T , equivalent to an instantaneous remineralization of terrigenous DOC at the river mouth, allows only 7% of that carbon to leave the Arctic (in the form of C T ), a smaller amount than estimated by Hansell et al. (2004), especially because their estimate does not consider the outflow of C T from remineralized terrigenous DOC. That comparison further emphasizes that our model overestimates the loss of terrigenous DOC to the atmosphere and underestimates the export of that dissolved carbon out of the Arctic. Specifying a semilabile pool of terrigenous DOC in the model would lead to greater offshore transport of DOC before it could be remineralized to C T . Consequently surface ocean pCO 2 and hence CO 2 outgassing would decrease near the river mouths and increase pCO 2 elsewhere in the Arctic Ocean as this semilabile terrigenous DOC would be remineralized later.
On the other hand, the GN2-based estimate for R DOC used as a model boundary condition underestimates the data-based estimate by 32% (Table 1). If the model-imposed R DOC had been larger, the outgassing of CO 2 to the atmosphere from doubling R DOC would also have been larger. Therefore, the model's low R DOC boundary condition compensates to some extent its overestimate of the loss of terrigenous DOC to the atmosphere caused by the direct remineralization of R DOC .
Unlike the effect from a doubling of R DOC , the doubling of R C T allows 96% of the added C T to remain in the ocean (approximately three fourths remains in the Arctic and approximately one fourth is found in the Atlantic after 70 years of simulation). Out of the C T that remains in the Arctic Ocean, 71% ends up in the central Arctic Ocean, while the remainder is distributed between the CAA (10%), the Barents Sea (6%), the Laptev Sea (5%), the Kara Sea (4%), the East Siberian Sea (3%), and the Chukchi and Beaufort Seas (1% each). Only 4% of that additional C T is emitted to the atmosphere (Table 5) because R A T is assumed to increase exactly in step with R C T . Therefore, the R ⋆ (1)) ratio increases in RIC while it decreases in ROC.
For simplicity, the PISCES model assumes that R A T equals R C T . In contrast, in three out of the five major Arctic rivers that drain directly into the Arctic Ocean, the A T :C T ratio varies from 0.91 to 0.94 . In the remaining two major Arctic rivers, the Kolyma and Ob, that ratio is 0.7 and 0.84, respectively. Thus, in the RIC simulation, the amount of added C T that is not buffered by an equivalent increase in A T , is up to 30%, a fraction that behaves as does the pure C T added in the ROC simulation where terrigenous DOC is assumed to be instantaneously remineralized to C T at the river mouth. Thus, by assuming R A T to be equal to R C T , we underestimate the outgassing of CO 2 in RIC and overestimate the corresponding increase in pH and Ω arag as detailed in the following section. (1)) from 0.71 to 0.55, i.e., for ΔROC, which in turn lowers the Arctic coastal ocean A T :C T ratio along with pH and Ω arag . Basinwide average changes of pH and Ω arag both reach −0.02, at most 7% of the respective changes due to the atmospheric CO 2 increase Global Biogeochemical Cycles 10.1029/2019GB006200 ( Table 6). Larger changes occur in the regional seas, for example, reaching up to −0.08 for pH and −0.09 for Ω arag , at most 30% of the respective changes due to the atmospheric CO 2 increase. The largest simulated changes occur very near to river mouths and are up to 3-37 times larger than surrounding regional averages ( Figure 5), as they decline sharply with distance away from each river mouth while the added C T (from R DOC ) is mostly lost via CO 2 outgassing (section 4.2). Thus, doubling R DOC mainly affects the waters very close to river mouths, where local CO 2 outgassing nearly reaches levels seen for basinwide average ingassing from the doubling of atmospheric CO 2 ( Figure 5). In the model, there is also some signature of low pH and low Ω arag along the transpolar drift ( Figures 5 and 6) from offshore transport of coastal waters into the central Arctic Ocean.

Ratio and Coastal Ocean Acidification
Contrary to case for R DOC , doubling the model's R C T comes along with an equal, simultaneous increase in R ⋆ A T by definition. Thus, doubling R C T increases the R ⋆ A T ∶ R ⋆ C T ratio from 0.55 to 0.71 because R DOC is nonzero (equation (1)) and in the model that is instantly converted to a flux of C T only. Furthermore, this increase in the R ⋆ A T ∶ R ⋆ C T ratio leads to an increase in the ocean A T :C T ratio. Hence, there are increases in the Arctic Ocean's surface average pH (0.02) and Ω arag (0.06) ( Table 6). The largest simulated regional average pH increases from doubling of R C T occur in the Laptev Sea (0.06) and Beaufort Sea (0.04). There that doubling cancels 21% and 15% of the corresponding pH declines from the atmospheric CO 2 doubling, based on our idealized assumption that atmospheric CO 2 , R C T , and R DOC all increase at the same rate. The largest changes in surface pH occur near river mouths and are 2 to 3 times larger than surrounding regional changes, a smaller factor than for the case of R DOC doubling where river mouth maxima were 3-37 times larger than corresponding regional changes. These lower river mouth coastal sea gradients from the doubling of R C T Global Biogeochemical Cycles 10.1029/2019GB006200 are attributed to the lack of CO 2 outgassing from increasing R C T and the transport of the additional C T and A T from rivers away from river mouths.
Regarding simulated Ω arag , its largest regional changes from the doubling of R C T occur in the Laptev Sea (0.19) and Beaufort Sea (0.11; Table 6), both regions with relatively higher R C T to R DOC ratios (Table 1). Those increases in those two seas offset 35% and 19% of the declines from doubling of atmospheric CO 2 . As seen for pH, the maximal changes of Ω arag occur very near river mouths and are 2 to 3 times larger than surrounding regional changes.
However, our results are subject to two simplifications concerning riverine fluxes: the instantaneous conversion of terrigenous DOC to C T and R A T :R C T =1. These two simplifications partly compensate one another, as the comparison between our R ⋆ A T :R ⋆ C T ratio and a reconstructed R ⋆ A T :R ⋆ C T ratio based on more realistic cases (R DOC lability of 50%, and R A T :R C T as in  suggest ( Figure S1). The two simplifications also cause the temporal change of the R ⋆ A T :R ⋆ C T ratio in the ROC and RIC simulations to be 26-55% larger than the more realistic cases in four out of five major Arctic rivers. Only for the Ob River do the two simplifications lead to a reduction (by 7%) of the change in R ⋆ A T :R ⋆ C T from the doubling of R DOC and R C T . Basinwide then, our estimates of the effects of doubling R DOC and R C T on pH and Ω arag are likely overestimated.

Sensitivity Factors
The sensitivity factors evolve over time but tend to stabilize by year 70 (Figures S2-S5). Sensitivity factors for carbon storage ( ) take longer to stabilize because the timescale for penetration of anthropogenic carbon into the deep ocean is much longer than our 70-year simulation, unlike for the shorter timescales of the other sensitivity factors, which are tied to the surface.
These sensitivity factors for air-to-sea CO 2 flux ( ), ocean carbon storage ( ), NPP ( ), and surface pH ( ) are summarized as time-integrated values at year 70 when forcing variables are doubled. Our simulated sensitivity of global ocean anthropogenic carbon uptake to increasing atmospheric CO 2 ( global ΔCO 2 ) may be compared to estimates from an ensemble of CMIP5 models forced under a 1%/year atmospheric CO 2 increase scenario (Arora et al., 2013). Our global ΔCO 2 of 0.9 Pg C/ppm lies at the upper limit of the CMIP5 model range (0.7-0.9 Pg C/ppm). The model sensitivity for the Arctic Ocean alone ( ΔCO 2 ) is 4.8 Tg C/ppm, only 0.5% of the analogous global sensitivity global ΔCO 2 even though the Arctic Ocean surface comprises 4% of the global ocean surface area. The relatively low ΔCO 2 for the Arctic Ocean is consistent with the regional distribution of the CMIP5 results (Ciais et al., 2013, Figure 6.22), which show a lower air-to-sea flux sensitivity to increasing atmospheric CO 2 for the Arctic Ocean compared to global values. In contrast, the sensitivity for Arctic Ocean carbon storage ΔCO 2 (22.3 Tg C/ppm) is 2% of the corresponding global sensitivity, although the Arctic Ocean contains only 1% of the global ocean's water volume. That 2:1 ratio is consistent with that seen for the data-based estimates of anthropogenic carbon storage in the Arctic Ocean (Tanhua et al., 2009). The difference between the Arctic ΔCO 2 and ΔCO 2 confirms that most of today's anthropogenic carbon stored in the Arctic Ocean is imported from the Atlantic and Pacific (Terhaar et al., 2019a).
In a consistent fashion, we calculated the sensitivity of the air-to-sea CO 2 flux and carbon storage to changes in riverine carbon input. The sensitivity of the air-to-sea CO 2 flux to changes in R DOC ( ΔROC ) is −0.86 Tg C (Tg C) −1 , indicating that 86% of additional C T from R DOC has been lost by CO 2 outgassing (Table 7). Conversely with the increase in R C T , only 5% of the additional C T is lost to the atmosphere ( ΔRIC = −0.05 Tg C (Tg C) −1 ). In terms of Arctic Ocean carbon storage, increasing R DOC enhances that by 0.05 Tg C for every teragram of carbon of R DOC that is added to the ocean as C T . This low regional ΔROC is in line with the corresponding largely negative ΔROC , which is consistent with most of the added C T from the increase in R DOC being lost to the atmosphere. Conversely, the Arctic's ΔRIC of 0.71 Tg C (Tg C) −1 is 14 times larger than its ΔROC (Table 7), highlighting that most of the additional C T from R C T remains in the ocean since there is an equal increase in R ⋆ A T and hence surface ocean pCO 2 is hardly affected (section 4.2). For both ΔROC and ΔRIC, the sum of the absolute numbers of and is below 1, indicating that after 70 years, the remainder Global Biogeochemical Cycles 10.1029/2019GB006200 (difference relative to 1.0) has left the Arctic Ocean laterally to the North Atlantic, that is, 0.09 Tg C (Tg C) −1 for ΔROC and 0.24 Tg C (Tg C) −1 for ΔRIC.
Likewise, we are interested in the sensitivities of NPP, the air-to-sea CO 2 flux, and ocean carbon storage to riverine nutrient input. The increase in Arctic Ocean NPP from the increase in the riverine nitrogen flux is ΔRUT = 1.07 Tg N (Tg N) −1 (Table 7), that is, more than is fueled from Arctic rivers alone. That extra Arctic NPP may be fueled in part by additional river nutrients added from outside the Arctic, given that R NUT is increased in all rivers globally and assuming there is sufficient transport from adjacent regions. It could also come in part from the remineralization of the additional organic matter produced by the increases in R NUT , for example, on the Arctic shelf. The increase in NPP also drives increases in the Arctic Ocean air-to-sea CO 2 flux ( ΔRUT = 2.29 Tg C (Tg N) −1 ) and carbon storage ( ΔRUT = 0.58 Tg C (Tg N) −1 ). The difference between ΔRUT and ΔRUT can be explained by lateral export to the Atlantic Ocean and the enhanced storage of carbon in organic matter and POC by enhanced NPP. Based on the model's C:N molar ratio of 122:16, an additional carbon uptake of 6.99 Tg C (Tg N) −1 would have been expected if every mole of carbon consumed by NPP would have been replaced via the air-to-sea CO 2 flux. Yet, the simulated ΔRUT reveals that only 33% of the C T removed by NPP is replaced via invasion of atmospheric CO 2 . This result is consistent with an earlier estimate of air-to-sea CO 2 flux enhancement from NPP by Orr and Sarmiento (1992) for the global ocean (43%), although here the focus is on the Arctic Ocean with its large shelf seas and high riverine inputs.
Lastly, sensitivities of surface ocean pH were assessed with respect to increases in atmospheric CO 2 , R DOC , and R C T . The sensitivity of Arctic surface pH to the increase in R DOC ( ΔROC ) is −0.029 (Pg C) −1 . Increasing R C T along with R ⋆ A T leads to an increase in pH by 0.011 (Pg C) −1 ( ΔRIC ) because the simultaneous and equal increases in R C T and R ⋆ (1)). Relative to these basinwide average sensitivities, locally sensitivities can be several times larger, for example, near river mouths (sections 4.1-4.3).
When integrated over the 70-year simulation, carbon storage from the doubling of riverine R DIC fluxes amounts to 18% of carbon storage from the doubling of atmospheric CO 2 (Table 7). In comparison, carbon storage from the doubling of riverine R DOC and R NUT fluxes each amount to 0.5% of carbon storage from the doubling of atmospheric CO 2 . In terms of the air-to-sea CO 2 flux, the doubling of R DOC causes an outgassing that amounts to 40% of the magnitude of the ingassing from the doubling of atmospheric CO 2 , while the doubling of R DIC has the same opposing effect but amounts to only 5%. Conversely, the doubling of R NUT increased the simulated air-to-sea CO 2 flux by an amount that is equivalent to 12% of that from the doubling of atmospheric CO 2 . Changes in mean surface pH from the doubling of riverine R DOC and R DIC fluxes each amount to 7% of the change from the doubling of atmospheric CO 2 .

Conclusions
This study offers a preliminary assessment of the extent to which certain key biogeochemical characteristics of the Arctic Ocean are sensitive to changes in river fluxes of carbon and nutrients. It provides a quantitative view that helps to disentangle how these characteristics are affected across the Arctic Ocean by riverine inputs of terrigenous organic carbon, inorganic carbon, and nutrients. Doubling the riverine nutrient flux increases Arctic NPP by 11% on average, by more than 30% in the Laptev and Beaufort Seas, and by up to 100% near river mouths. Doubling riverine DOC fluxes enhances ocean CO 2 outgassing, nearly offsetting the influx of anthropogenic CO 2 through the Arctic Ocean's air-sea interface from an imposed simultaneous doubling of atmospheric CO 2 . However, much more anthropogenic carbon enters the Arctic Ocean laterally. Doubling riverine DOC fluxes also reduces the R ⋆ A T :R ⋆ C T ratio, thus enhancing ocean acidification. Average changes in surface pH and Ω arag due to the doubling of riverine DOC fluxes amount to at most 7% of the changes due to the doubling of atmospheric CO 2 , while that proportion reaches up to 30% in regional seas and up to 100% close to river mouths. Conversely, doubling river C T fluxes increases the R ⋆ A T : R ⋆ C T ratio and thus reduces ocean acidification. That reduction is at most 8% of the opposite effect from the doubling of atmospheric CO 2 in terms of the basinwide average, but it amounts to up to 21% of regional sea averages and up to 50% near river mouths.
This study takes a first step toward assessing the influence of changing riverine input on Arctic Ocean biogeochemistry. The approach is limited by model simplifications, for example, its relatively coarse lateral resolution (1 • ), which does not explicitly resolve mesoscale and submesoscale processes (Holt et al., 2008). Beyond eddies, coastal upwelling and internal tides are poorly represented (Holt et al., 2014;Nurser & Bacon, 10.1029/2019GB006200 2014. Such weaknesses also affect simulated ocean biogeochemistry (Holt et al., 2008), for example, the air-to-sea CO 2 flux (Bianchi et al., 2005). Future work with nested models could potentially alleviate these physical concerns while limiting the computational investment.
These results suggest that the effects of riverine DOC fluxes should also be considered in the debate about whether the Arctic Ocean will become a source or a sink of CO 2 (Bates & Mathis, 2009;Cai et al., 2010;Manizza et al., 2019;Roy et al., 2011). Likewise, riverine nutrient fluxes should be accounted for in another debate, the one about how nutrient supply to the surface Arctic Ocean will change and whether those changes will enhance or reduce NPP in the future (Arrigo & van Dijken, 2015;Cai et al., 2010;Vancoppenolle et al., 2013). Indeed, the wide divergence between NPP projections among the CMIP5 models (Vancoppenolle et al., 2013) may be partly due to the lack of accounting of riverine nutrient fluxes in some models. Not accounting for the effects of riverine carbon fluxes may also lead to biased projections of pH and Ω arag in Arctic shelf seas. Our finding that doubling riverine C T fluxes reduces coastal ocean acidification is counterintuitive if one does not consider the implicit simultaneous increase in A T . Although this finding is based on two simplifications, more realistic assumptions still lead to a future increase in R ⋆ A T :R ⋆ C T and hence reduced acidification from an increase in riverine inorganic carbon delivery.
Overall, our results suggest that even in the Arctic Ocean, where riverine inputs are proportionally the largest, the effects of decadal-to-centennial changes in river fluxes of carbon and nutrients on open ocean biogeochemistry are relatively small. Yet they also suggest that for coastal seas, the influence of riverine input of carbon and nutrients should not be neglected, as often is the case in ocean models.