Oceanic CO2 outgassing and biological production hotspots induced by pre-industrial river loads of nutrients and carbon in a global modelling approach

Rivers are a major source of nutrients, carbon and alkalinity to the global ocean. In this study, we firstly estimate pre-industrial riverine loads of nutrients, carbon and alkalinity based on a hierarchy of weathering and terrestrial organic matter export models, while identifying regional hotspots of the land-ocean exports. Secondly, we implement the riverine loads into a global biogeochemical ocean model and describe their implications for oceanic nutrient concentrations, the net primary production (NPP) and air-sea CO2 fluxes globally, as well as in a regional shelf analysis. Thirdly, we quantitatively assess 5 the terrestrial origins and the long-term oceanic fate of riverine carbon in the framework. We quantify annual pre-industrial riverine loads of 3.7 Tg P, 27 Tg N, 158 Tg Si and 603 Tg C. We thereby identify the tropical Atlantic catchments (20% of global C), Arctic rivers (9% of global C) and Southeast Asian rivers (15% of global C) as dominant providers of carbon to the ocean. The riverine exports lead to a net global oceanic CO2 source of 231 Tg C yr-1 to the atmosphere, which mainly results from inorganic riverine carbon loads (183 Tg C yr-1), and from organic riverine carbon inputs (128 Tg C yr-1). Additionally, 10 a sink of 80 Tg C yr-1 is caused by the enhancement of the biological carbon uptake by dissolved inorganic nutrient inputs and the resulting alkalinity production. While large outgassing fluxes are mostly simulated in proximity to major river mouths, substantial outgassing fluxes can also be found further offshore, most prominently in the tropical Atlantic. Furthermore, we find evidence for the interhemispheric transfer of carbon in the model; we detect a stronger relative outgassing flux (49% of global river induced outgassing) in the southern hemisphere in comparison to the hemisphere’s relative riverine inputs 15 (33% of global river inputs), as well as an outgassing flux of 17 Tg C yr-1 in the Southern Ocean. Riverine loads lead to a strong increase in NPP in the tropical West Atlantic, Bay of Bengal and the East China Sea (+166%, +377% and +71% respectively). While the NPP is not strongly sensitive to riverine loads on the light-limited Arctic shelves, the CO2 flux is strongly altered due to substantial dissolved carbon supplies to the region. While our study confirms that the ocean circulation is the main driver for open ocean biogeochemical distributions, it reveals the necessity to consider riverine exports for the 20 representation of heterogeneous features of the coastal ocean and to represent riverine-induced carbon outgassing in the ocean. It also underlines the need to consider the long-term volcanic CO2 flux to close the atmospheric carbon budget in a coupled land-ocean-atmosphere setting.


Introduction
Rivers deliver substantial amounts of carbon (C), phosphorus (P), nitrogen (N), silica (Si), iron (Fe) and alkalinity (Alk) to the ocean (Seitzinger et al., , 2010Dürr et al., 2011;Beusen et al., 2009;Tréguer and De La Rocha, 2013;Beusen et al., 2016). In the ocean, these compounds undergo biogeochemical transformations and are ultimately either buried in the sediment or are outgassed to the atmosphere (Froelich, 1988;Sarmiento and Sundquist, 1992;Aumont et al., 2001;Stepanauskas 5 et al., 2002;Dagg et al., 2004;Krumins et al., 2013). In global ocean biogeochemistry models, riverine inputs and their contributions to the oceanic cycling of C have been strongly simplified or even completely omitted. In this study, we attempt to improve the understanding of the long-term effects of riverine loads in the ocean by firstly estimating the magnitudes of biogeochemical riverine loads for the pre-industrial time period and secondly by assessing their long-term implications for marine biogeochemical cycles in a global ocean model. 10 Natural riverine C and nutrients originate from the weathering of the lithosphere and from organic matter exports of the terrestrial biosphere (Ludwig et al., 1998). Weathering directly releases nutrients (P, Si and Fe) that can be taken up by the terrestrial ecosystems, or exported directly to aquatic systems . In these ecosystems, these nutrients are reported to enhance the carbon uptake due to their limitation of the biological primary production (Elser et al., 2007;Fernández-Martínez et al., 2014). During the weathering process, Alk and dissolved inorganic C (DIC) are released, while 15 CO 2 is drawn down from the atmosphere (Amiotte Suchet and Probst, 1995;Meybeck and Vörösmarty, 1999;Hartmann et al., 2009). Spatially explicit weathering models could potentially be used to provide the weathering release yields of nutrients, Alk and C (e.g. Hartmann et al. (2014)), as well as to quantify the weathering-induced drawdown of atmospheric CO 2 in Earth System Models (ESMs, Ludwig et al. (1998); Hartmann et al. (2009); Roelandt et al. (2010); Hartmann et al. (2014)). These approaches rely on estimating chemical weathering rates as a first-order function of hydrology, lithology, rates of physical 20 erosion, soil properties and temperature (Amiotte Suchet and Probst, 1995;Hartmann et al., 2009Hartmann et al., , 2014. The terrestrial biosphere also provides C and nutrients (P, N, Si and Fe) to rivers. This happens dominantly through organic matter exports (Meybeck and Vörösmarty, 1999;Seitzinger et al., 2010;. The formation of organic matter through biological primary production on land consumes atmospheric CO 2 (Ludwig et al., 1998). Then, dissolved and particulate organic matter can be mobilized in soils and peatlands through leaching and physical erosion, and exported to freshwaters. 25 While the natural P and Fe within the terrestrial organic matter originate from weathering, C and N mostly originate from atmospheric fixation (Meybeck and Vörösmarty, 1999;Green et al., 2004).
Dissolved inorganic nutrients enhance the biological primary production in the ocean, which results in a net uptake of atmospheric CO 2 (Tyrrell, 1999). In turn, riverine inputs of DIC result in C outgassing in the ocean (e.g. Hartmann et al. (2009)). The remineralization of terrestrial organic matter provided by rivers releases DIC and nutrients to the ocean. The 30 composition of this organic matter and the time-scales of its degradation in the ocean have however been open questions for over two decades (Ittekkot, 1988;Hedges et al., 1997;Cai, 2011;Lalonde et al., 2014). In the case of terrestrial dissolved organic matter (tDOM), its reported degradation during its transit in rivers induce high C to nutrients ratios of its loadings to the ocean (i.e. C:P weight ratios of over 500, Meybeck (1982); Seitzinger et al. (2010)), which is likely due to the preference of organic matter rich in nutrients for bacterial metabolic activities (Hedges et al., 1997)). The degradation of tDOM is suggested to cause C outgassing in the coastal ocean due to its large release of C and low release of nutrients to the water column (Cai, 2011;Müller et al., 2016;Aarnos et al., 2018). Although the previous degradation of tDOM in rivers implies low biological reactivity of tDOM in the ocean, tDOM is not a major constituent of organic mixtures in open ocean seawater or in sediment pore water (Ittekkot, 1988;Hedges et al., 1997;Benner et al., 2005). It is thus reactive to a certain degree in the ocean, with 5 recent studies suggesting the abiotic breakdown of tDOM as a possible explanatory mechanism (Fichot and Benner, 2014;Aarnos et al., 2018). In the case of riverine particulate organic matter (POM), stronger gaps of knowledge exist (Cai, 2011).
Riverine POM has however been reported to affect the coastal ocean biogeochemistry regionally by controlling the availability of nutrients through its remineralization (Froelich, 1988;Dagg et al., 2004;Stramski et al., 2004). Furthermore, a substantial proportion of weathered P is exported to the ocean bound to iron-manganese oxide and hydroxide particle surfaces (Fe-P, 10 Compton et al. (2000)). Within the ocean, a fraction of P in Fe-P is thought to be desorbed, and thus converted to bioavailable dissolved inorganic phosphorus (DIP).
Pre-industrial P and N riverine loads likely strongly differed to their present-day loads due to dramatic anthropogenic perturbations of inputs to freshwaters during the 20 th century (Seitzinger et al., 2010;Beusen et al., 2016). Riverine exports of P and N are also suggested to have already been perturbed prior to 1850 due to increased soil erosion from land-use changes, to 15 fertilizer use in agriculture and to sewage sources (Mackenzie et al., 2002;Filippelli, 2008;Beusen et al., 2016). Since global modelling studies of the climate are usually initialized for 1850 (Giorgetta et al., 2013) or 1900 (Bourgeois et al., 2016), these exports from non-natural sources should also be taken into account in initial pre-industrial model simulations. While anthropogenic perturbations of C and Si have also been reported in published literature, they are less substantial at the global scale than for P and N (Seitzinger et al., 2010;Maavara et al., 2014Maavara et al., , 2017. 20 Until now, riverine point sources of biogeochemical compounds have been omitted or poorly represented in global ocean biogeochemical models, despite being suggested to strongly impact the biogeochemistry of coastal regions (Froelich, 1988;Stepanauskas et al., 2002;Dagg et al., 2004) and to cause a pre-industrial source of atmospheric CO 2 in the ocean (Sarmiento and Sundquist, 1992;Aumont et al., 2001;Gruber et al., 2009;Resplandy et al., 2018). This background CO 2 outgassing flux of 0.2 to 0.8 Gt C yr -1 is significant in the context of the present-day oceanic C uptake of around 2.3 Gt C yr -1 (IPCC, 2013). In 25 a modelling study, Aumont et al. (2001) derived terrestrial C loads from an erosion model and analyzed the oceanic outgassing caused by these C inputs. The impacts of riverine nutrients and Alk were however not considered and C was only added to the ocean as DIC, thus omitting terrestrial organic matter dynamics in the ocean. Da Cunha et al. (2007) analyzed the impacts of present-day riverine loads on the oceanic NPP in an ocean biogeochemical model covering the Atlantic Ocean for an analysis period of 10 years, a time period that likely does not suffice to assess their long-term implications on the biogeochemistry of 30 basin areas remote to land. Bernard et al. (2011) added present-day biogeochemical riverine loads to an global ocean model to focus on their implications for global opal export distributions, ignoring other impacts of the riverine inputs in the model. In a global ocean model as well, Bourgeois et al. (2016) quantified the coastal anthropogenic CO 2 uptake, but the impacts caused by biogeochemical riverine loads were not represented. To our knowledge, a study has yet to give a comprehensive overview of the magnitudes of pre-industrial riverine exports and their long-term implications for the global oceanic biogeochemical cycling in 35 a 3-dimensional model. An initial ocean state that considers pre-industrial riverine supplies has not been achieved in a global ocean model, which would be needed to investigate the temporal impacts associated with their perturbations. Past ocean model studies have also not represented the composition and oceanic dynamics of terrestrial organic matter, which differs to oceanic organic matter. Furthermore, it is often unclear under which criteria the Alk supplies to the ocean were constrained in many global ocean models. Regional coastal responses to the addition biogeochemical riverine loads have not been assessed at the 5 global scale until now, largely due to the incapability of the global ocean models to represent plausible continental shelf areas in the past (Bernard et al., 2011).
To address the knowledge gaps listed above, we 1. implement a representation of pre-industrial riverine loads into a global ocean biogeochemical model, considering both weathering and non-weathering sources of nutrients, C and Alk. We compare our estimates with a wide range of published literature values, while also determining regions of disproportionate contributions 10 to global terrestrial exports. 2. The implications of riverine inputs for oceanic net primary production (NPP) and air-sea CO 2 fluxes are assessed globally, as well as regionally in an analysis of shallow shelf regions. 3. We evaluate the origins and fate of riverine C quantitatively, while assessing the balance between the land C uptake and the oceanic C outgassing in the individual models. This balance of the land C uptake and its oceanic outgassing is then used to assess the potential implementation of riverine fluxes in a fully coupled land-atmosphere-ocean setting. 15

Methods
To address the objectives of this study, we derived the most relevant pre-industrial (1850) riverine loads of biogeochemical compounds to the ocean in dependence of pre-industrial ESM variables. The derived riverine loads were then incorporated into a global ocean biogeochemical model in order to assess their global and regional impacts in the ocean. In order to quantify the response of coastal regions to the riverine supplies, we also defined 10 regions with ocean depths of less than 250m for 20 analysis.

Deriving pre-industrial riverine loads
We focused on the riverine exports of P, N, Si, Fe, C and Alk at the global scale. The river catchments were defined by using the largest 2000 catchments of the Hydrological Discharge (HD) model (Hagemann and Dümenil, 1997;Hagemann and Gates, 2003), a component of the Max Planck Institute Earth System Model (MPI-ESM). They were derived from the runoff flow 25 directions of the model at a horizontal resolution of 0.5 degrees. The exorheic river catchments (catchments, which discharge into the ocean) were considered to be catchments which had river mouths at a distance of less than 500 km to the coastline in the HD model. Catchments with larger distances were considered to be endorheic catchments (catchments, which do not discharge into the ocean). The biogeochemical compound yields in these catchments were assumed to be retained permanently, whereas the riverine exports of exorheic catchments were added to the ocean. 30 We considered both weathering as well as non-weathering sources of P, N, Si, Fe, C and Alk to river catchments (Figure 1). These were derived, if possible, from spatially explicit models. Within the catchments, we accounted for transformations of P, N and Fe to organic matter through biological productivity on land and in rivers. These transformations were derived from globally fixed ratios with respect to organic C. The organic C in tDOM and POM was, in turn, assumed to ultimately originate from terrestrial biological CO 2 uptake (Ludwig et al., 1998) and was therefore not subtracted from the catchment DIC pools.
Additionally, a fraction of the catchment P was also assumed to have been adsorbed to Fe-P. Everything considered, rivers deliver terrestrial organic matter (tDOM and POM), inorganic compounds (DIP, Fe-P, DIN, DSi, DFe and DIC) and Alk to the 5 ocean in our approach (Figure 1).
The surface runoff, surface temperature and precipitation data used to drive the weathering release and land export models were obtained from output of a coupled MPI-ESM pre-industrial simulation (Giorgetta et al., 2013). We thereby used the annual 100 year means of pre-industrial data computed at a horizontal grid resolution of 1.875 degree. The runoff data was scaled by a factor of 1.59 to account for the runoff model bias with regard to global estimations (Fekete et al., 2002), a factor 10 that is discussed in Appendix B.

Terrestrial dissolved and particulate organic matter characteristics
We assumed that the pre-industrial loads of tDOM and POM did not strongly differ to their present-day loads at the global scale. Seitzinger et al. (2010), who modelled anthropogenic perturbations to riverine loads for 1970 and 2000, only simulated small changes in global POC and DOC loads over this time period.  suggest a total anthropogenic perturbation 15 of the global riverine C export to the ocean of around 10%, for which we did not account in this study.
The riverine loads of tDOM and POM were therefore directly derived from the DOC and POC river loads for the reference year 1970 (NEWS2), assuming no change between this time segment and the pre-industrial time period. These organic loads were determined from the models of Harrison et al. (2005) and Beusen et al. (2005). The Harrison et al. (2005) model quantifies the DOC catchment yields as a function of runoff, wetland area and consumptive water use. Beusen et al. (2005) describe the 20 POC catchment yields as a function of catchment suspended solids yields, which depend on grassland and wetland areas, precipitation, slope and lithology.
The tDOM and POM exports were assumed to consist of globally constant fractions of organic C, P, N and Fe. The tDOM C:P ratio was based on a C:P mole ratio of 2583:1 derived from Meybeck (1982) and Compton et al. (2000). The tDOM total N:P mole ratio was chosen to be 16:1, in order to represent the same N:P ratio as of the organic matter export to the sediment in 25 the ocean biogeochemistry model, which was also the reasoning for choosing the P:Fe mole ratio of 1:3.0 10 -4 . The C:N:P:Fe mole ratio of tDOM was therefore 2583:16:1:3.0 10 -4 . The C:P ratio of riverine POM is highly uncertain, but global C:P mole ratios from observational data (56-499) (Meybeck, 1982;Ramirez and Rose, 1992;Compton et al., 2000) suggest a much closer ratio to the production and export ratios of oceanic organic matter (mole C:P = 122:1, Takahashi et al. (1985)) than for tDOM. Due to this and current gaps of knowledge on the composition of riverine POM, we chose the same C:N:P:Fe ratio as

6
We derived the P weathering yields from a spatio-temporal model , which quantifies the P weathering release in relation to the SiO 2 and cation release induced by weathering. The model was originally calibrated for the extensive data set of 381 river catchments of the Japanese Archipelago with runoff and lithological characteristics (Hartmann and Moosdorf, 2011) being used as sole model variables. The model was then corrected for temperature and soil shielding effects at the global scale : where F Prelease is the chemical weathering rate of P per area (t km -2 yr -1 ), b P,i is an empirical factor representing the rate of P weathering of the lithology i, q is the surface runoff (mm yr -1 ), F i (T) is a lithology-dependent temperature function and F S is a parameter to account for soil shielding.
The lithology types consisted of 16 lithological classes from the global lithological database GliM  dorf, 2012). The lithology-dependent factor b P,i represents the chemical weathering rate factor for SiO 2 + cations (b SiO2+Cat ) multiplied with the relative P content (b Prel,i ): The parameters b SiO2+Cat,i and b Prel,i for each lithology i can be found in Hartmann et al. (2014).
The temperature correction function F(T) is an Arrhenius relationship for basic (rich in iron and magnesium) and acid (high 15 silica content) lithological classes, with activation energies normalized to the average temperature of the calibration catchments of the study . For acid rock lithologies, an activation energy of 60 kJ/mole was assumed, whereas for basic rock types 50 kJ/mole was used. Carbonate lithologies do not have a temperature correction function due to the absence of a clear relationship to field data, as well as due to uncertainties in the mechanisms of a temperature effect on carbonate weathering Romero-Mujalli et al., 2018). 20 A soil shielding factor F S was considered to represent the inhibition of weathering by certain types of soils. These soils, which have low physical erosion rates, develop a chemically depleted thick layer, shielding them from water supply and thus partly preventing the weathering of the soil aggregates (Stallard, 1995). An average soil shielding factor of 0.1 was estimated for the soils Ferrasols, Acrisols, Nitisols, Lixisols, Histosols and Gleysols while performing a calibration at the global scale in Hartmann et al. (2014). 25

Non-weathering P sources
Since the P cycle was already perturbed in the assumed pre-industrial state (1850) due to anthropogenic activities (Mackenzie et al., 2002;Filippelli, 2008;Beusen et al., 2016), we also accounted for P sources other than weathering (P nw,catch ). Similarly to Beusen et al. (2016), we derived the non-weathering catchment source of P as the sum of fertilizer (P fert,catch ), sewage (P sew,catch ) and allochthnonous P inputs (P alloch,catch ):

30
P nw,catch = P fert,catch + P sew,catch + P alloch,catch P fert,catch was the P input to river catchments originating from the agricultural application of fertilizers, manure and organic matter (1.6 Tg P yr -1 globally, from Beusen et al. (2016)), P sew,catch was the P input from sewage (0.1 Tg P yr -1 globally, Morée et al. (2013)) and P alloch,catch represented allochthonous organic matter P inputs (1 Tg P yr -1 , which were simplified as vegetation in floodplains in Beusen et al. (2016)). All of these P inputs were estimated for year 1900 due to previous estimates not being available to our knowledge. Since our framework was developed with the aim of being used in ESMs, we 5 assumed soil equilibrium since this is the initial state criteria used in state-of-the-art model simulations. Therefore, P exports due to perturbations of organic matter erosion in soils were not considered. The global distribution of anthropogenic P inputs (agricultural+sewage) to catchments was assumed to be the same as the distribution of contemporary anthropogenic P inputs, which was derived from the NEWS2 study: For each catchment, we estimated the annual P exports of the individual catchments (P total,catch ), as the sum of the catchment weathering yields (P w,catch ) and of non-weathering catchment P inputs (P nw,catch ): P total,catch = P w,catch + P nw,catch P total,catch was then fractionated into inorganic P (IP catch ) and organic P bound within tDOM (DOP) and POM (POP), which was assumed to have been taken up on land or in freshwaters by the biology at the prescribed fixed organic C:P ratio. The 20 remaining catchment P was assumed to be IP catch : The IP was thereafter fractionated into DIP and Fe-P (IP bound to iron-manganese oxide and hydroxide) at a ratio r inorg (DIP:Fe-P = 1:3, approximated from the pre-industrial global export estimates of Compton et al. (2000)) and (1-r inorg ): 25 and: We did not consider detrital particulate inorganic P exports from the physical weathering of rock. This P is not bioavailable in rivers or in the coastal ocean due to strong covalent bonds of the detrital mineral structures (Compton et al., 2000).

Nitrogen and iron
The N inputs to river catchments were derived from a globally fixed mole N:P ratio of 16:1 (Takahashi et al., 1985) with regard to P inputs for all species (DIN:IP, N:P in tDOM and N:P in POM). In published literature, higher N:P ratios for all riverine species exports are suggested, with Beusen et al. (2016) simulating a pre-industrial N:P mole ratio of 21:1 and Turner et al. (2003) reporting higher N:P ratios of over 16:1 in a synthesis of measurements in major rivers. However, denitrification in 5 river estuaries, which is not taken into account in the mentioned studies, also removes N (3-10 Tg N yr -1 , Seitzinger et al. (2005)). According to our calculations, this estuarine N sink could reduce the N:P ratio of the global exports to around 16:1 (Supplementary Information S.1.3). For Fe, we used a P:Fe mole ratio of 1:3.0 10 -4 , which is the Fe:P export ratio of organic matter in the ocean biogeochemistry model, to quantify Fe inputs to catchments for all species (IP:DFe, P:Fe in tDOM and P:Fe in POM). The iron of the Fe-P iron-manganese oxides and hydroxides was not considered to be bioavailable.

Dissolved inorganic carbon and alkalinity
The DIC and Alk weathering release, as well as the draw down of atmospheric CO 2 induced by weathering, were quantified according to the studies of Hartmann et al. (2009) and Goll et al. (2014). Weathering reactions cause an uptake of atmospheric CO 2 and a release of DIC as HCO 3 during carbonate weathering: 15 and silicate weathering: The equations (9) and (10) dictate the release of 1 HCO 3 -(thus 1 DIC and 1 Alk) for each mole of CO 2 taken up during the weathering of silicate lithologies, and of 2 HCO 3 -(thus 2 DIC and 2 Alk) for the uptake of each mole of CO 2 drawn down during the weathering of carbonate lithologies. 20 The release equations from Hartmann et al. (2009) and Goll et al. (2014) quantify the lithology (i) dependent HCO 3 weathering release as a function of surface runoff (q), temperature (F i (T)), soil shielding F S and a weathering parameter b C,i : The parameter b C,i is dependent on the weathering rate of the lithology and the composition of the lithology.
The catchment DIC and Alk catchment loads were the HCO 3 weathered annually within the catchments, assuming con- 25 servation of Alk along the land-ocean continuum (Amiotte Suchet and Probst, 1995). We therefore assumed a DIC:Alk loads ratio of 1:1, without considering additional DIC sources from respiration of organic matter in soil pore water, groundwater or in rivers. River observational data however show that the riverine HCO 3 and total DIC mole exports rarely deviate by more than 10% (Meybeck and Vörösmarty, 1999;Araujo et al., 2014).

Silica
To quantify the global spatial distribution of riverine DSi exports, we used the model developed by Beusen et al. (2009): where F DSiO2 is the export of DSi in Tg SiO 2 yr -1 km -2 , ln(prec) is the natural logarithm of the precipitation in mm d -1 , volc is the area fraction covered by volcanic lithology (no dimension), bulk is the bulk density of the soil in Mg m -3 , slope is the 5 average slope in m km -1 , and b prec , b volc , b bulk and b slope are the regression coefficients estimated in Beusen et al. (2009). For the precipitation, we used the pre-industrial model output from the MPI-ESM, whereas the volcanic area originated from Dürr et al. (2005), the soil density from Batjes (1997) and Batjes (2002), and the average slope from the Global Agro-Ecological Zones database (FAO/IIASA). The exports were aggregated for the HD model catchment areas. The loads that were generated by applying the Beusen et al. (2009) model were also converted to Tg Si and are given accordingly in our study. We approx-10 imated the weathering release of Si using this DSi model output, therefore assuming that sinks of Si on land or in rivers are compensated by sources from these ecosystems (a further equilibrium assumption). We also neglected the riverine exports of particulate Si originating from physical erosion.  (2006)). The atmospheric forcing data set is a mean annual cycle for atmospheric parameters at a daily time step constructed from ECMWF reanalysis project ERA40 data. For the freshwater forcing climatology, the HD model was driven with the ERA40 atmospheric data (Hagemann and Gates, 2001). 25 The Hamburg Ocean Carbon Cycle model (HAMOCC) simulates the inorganic carbon chemistry, biological transformations, nutrient cycling and gas dynamics in the oceanic water column, sediment and at the air-sea interface. The model core along with its equations are described in Ilyina et al. (2013), but we also accounted for more recent modifications explained in Mauritsen et al. (2018). These included incorporating dynamic nitrogen fixation by cyanobacteria (Paulsen et al., 2017), following recommendations of the OMIP protocol (Orr et al., 2017) and including nitrogen deposition according to the Coupled The inorganic carbon chemistry is based on Maier-Reimer and Hasselmann (1987), with adjustments in the calculation of chemical constants as described in the OMIP protocol (Orr et al., 2017). Total DIC and total Alk are thereby prognostic tracers from which the carbonate species are determined diagnostically. Total DIC includes all carbonate species and total Alk includes 5 carbonate as well as borate alkalinity.

Ocean
The dynamics of organic matter cycling in the ocean are based on an extended NPZD type model with pools of nutrients, phytoplankton, zooplankton, detritus (POM), DOM and cyanobacteria. The phytoplankton growth follows Michaelis-Menten kinetics as a function of temperature, light and nutrient availability. A constant ratio (C:N:P:Fe = 122:16:1:3.0 10 -4 , Takahashi et al. (1985)) dictates the composition of all oceanic organic matter. Both oceanic DOM and POM are advected according to 10 the ocean physics, and the POM also sinks as a function of depth (Martin et al., 1987). The aerobic remineralization of organic matter takes place when the O 2 concentration is above a threshold, whereas at low enough O 2 concentrations, denitrification and sulfate reduction can take place. Furthermore, the phytoplankton produce opal shells when dissolved silica is available, and CaCO 3 shells when dissolved silica is depleted. The CaCO 3 and opal sink at constant rates.
HAMOCC also contains a 12 layer sediment module that simulates the same remineralization and dissolution processes 15 as in the water column for sediment constituents (Heinze et al., 1999). The sediment consists of a fraction of pore water, which contains dissolved inorganic compounds (e.g. DIC and DIP). POM, CaCO 3 and opal fluxes from the water column are deposited to the top sediment layer. There is a diffusive inorganic compound flux at the water-sediment interface and a particulate flux from the bottom layer to a diagenetically consolidated burial layer.
Dust is added through atmospheric deposition according to input fields of Mahowald et al. (2006). The model considers the 20 dynamics of CO 2 , O 2 , N 2 O and N 2 and their exchange at the ocean-atmosphere interface. Since we model a pre-industrial state in this study, we used constant atmospheric concentrations for CO 2 of 278 ppmV (Etheridge et al., 1996).

Treatment of the river loads in the ocean biogeochemistry model
The biogeochemical riverine loads were added to the ocean surface layer in HAMOCC constantly over the whole year. The locations of major river mouths were corrected manually on a case to case basis for large rivers in order to reproduce the same 25 locations as of the OMIP freshwater inputs.
The riverine dissolved inorganic compounds (DIC, Alk, DIP, DIN, DSi, DIC, DFe) were added to the HAMOCC in their species pools. We added 80% of P contained in the riverine Fe-P to the oceanic DIP pool, in order for the amount of bioavailable Fe-P to be comparable to the estimated range given in Compton et al. (2000) (1.1-1.5 Tg P yr -1 ). The rest of the Fe-P pool was considered to be unreactive in the ocean. The riverine POM was added to the oceanic POM pool in HAMOCC. 30 For tDOM, we extended HAMOCC with a new model tracer that was characterized by the described C:N:P:Fe mole ratio.
tDOM was mineralized as a function of the tDOM concentration at a rate constant k rem,tDOM and of an oxygen limitation factor (Γ O2 ), which decreases the maximum potential remineralization rate dependently on the O 2 concentration: Since a large fraction of tDOM delivered by rivers is already strongly degraded, it is to a certain extent resistant to microbial degradation in the ocean (Ittekkot, 1988;Vodacek et al., 2003). We therefore assumed a slightly lower remineralization rate constant of tDOM (k rem,tDOM ) compared to oceanic DOM (0.003 versus 0.008 d -1 for oceanic DOM). This rate constant is function used (Γ O2 ) was the same as the function used for oceanic DOM, which is described in Mauritsen et al. (2018).

Pre-industrial ocean biogeochemistry model simulations
We performed two ocean model simulations in order to assess the impacts of biogeochemical riverine fluxes in terms of their were approximately stable, we used the 100-year mean of these fluxes to perform a 10'000-year simulation of the sediment sub-model. The resulting sediment state was then coupled back to the water column for a simulation of 2'000 more years.
We used REF as a reference simulation in order to compare RIV to, which enabled us to compare the impacts of riverine loads at plausible pre-industrial magnitudes and locations (RIV) to REF, where the vast majority of inputs were added directly to the open ocean. For the analysis of the resulting ocean biogeochemistry, we used the last 100-year means of RIV and REF.

Definitions of coastal regions for analysis
To investigate the impacts of riverine exports on coastal regions, we selected 10 shallow shelf regions characterized by large 5 riverine loads that cover a variety of latitudes (Table 1). The shelves were defined to have depths of less than 250m. The cutoff sections perpendicular to the coast were done according to MARgins and CATchements Segmentation (MARCATS) (Laruelle et al., 2013), except for the Beaufort Sea, Laptev Sea, North Sea and Congo shelf, where we used the COastal Segmentation and related CATchments (COSCATs) definitions (Meybeck et al., 2006) due to the vastness of their MARCATS segmentations.  (2006)). The shelf classes were defined as in Laruelle et al. (2013).

Global weathering release
The weathering release models provide global annual means for the pre-industrial weathering release of P, Si, DIC and Alk (Table 2), as well as of their spatial distributions (Figure 2). The calculated global P release is 1.34 Tg P yr -1 when considering the runoff scaling factor of 1.59, which fits within the P release range of 0.8 -4 Tg P yr -1 reported in published literature (Compton et al., 2000;Wang et al., 2010;Goll et al., 2014;Hartmann et al., 2014). While the Hartmann et al. (2014) and 15 Goll et al. (2014) estimates (1.1 and 0.8-1.2 Tg P yr -1 , respectively) originate from the same weathering model as is used here, the 1.9 Tg P yr -1 value reported in Wang et al. (2010) was constructed by upscaling measurement data points. In a further study, Compton et al. (2000) provide a quantification of the prehuman phosphorus cycle while distinguishing between landocean fluxes from shale-erosion as well as from weathering. The reported total prehuman riverine P export with a weathering source is given by averaging non-detrital P species concentrations from unpolluted rivers and multiplying them with global 5 runoff estimates, thereby yielding an estimate of 2.5 -4.0 Tg P yr -1 . Both estimates that originate from the upscaling of river measurements are therefore higher than the P weathering fluxes provided in the modelling approach in this study, in Goll et al. (2014) and in Hartmann et al. (2014), which suggests further effort might be needed to constrain the global P weathering release more accurately.
For Si, we model a global release of 168 Tg Si yr -1 . This is within the range estimated by Beusen et al. (2009) (158-199 Tg 10 Si yr -1 ), who used the same model while using present-day observational data to drive the model for precipitation. Our estimate is also comparable with the 173 Tg Si yr -1 value reported by Dürr et al. (2011).
The modelled DIC and Alk release amounts to 374 Tg C yr -1 / 37.8 Teq yr -1 globally, which also takes into account the runoff scaling factor. By extrapolating from measurement data of 60 large river catchments, Meybeck (1982) suggests that the global DIC export to the ocean is around 380 Tg yr -1 and originates directly from weathering. Further modelling studies 15 also provide similar weathering release estimates of 260 to 300 Tg C yr -1 (Berner et al., 1983;Amiotte Suchet and Probst, 1995). The atmospheric CO 2 drawdown induced by weathering is directly related to the weathering release of DIC since silicate weathering draws down 1 mole of CO 2 per mole HCO 3 -, and carbonate weathering draws down 0.5 mol of CO 2 per mole HCO 3 released (Eq. (9) and Eq. (10)). While we model a CO 2 drawdown flux of 280 Tg C yr -1 induced by weathering, drawdown fluxes of 220-440 Tg C yr -1 are reported in published literature (Gaillardet et al., 1999;Amiotte Suchet et al., 2003;20 Hartmann et al., 2009;Goll et al., 2014). The results imply that, of the 374 Tg C yr -1 DIC released by weathering, 280 Tg C yr -1 originates from atmospheric C drawdown, while the rest originates from the weathering DIC release from the carbonate lithology (94 Tg C yr -1 ).
We observe hotspots that contribute disproportionately to the weathering release yields, with the Amazon, Southeast Asia as well as Northern Europe and Siberia acting as strong sources of weathering (Figure 2). The disproportionate importance 25 of these regions for the global weathering release is largely due to wet and/or warm climates, as well as to their lithological characteristics (Amiotte Suchet and Probst, 1995;Gaillardet et al., 1999;Hartmann and Moosdorf, 2011).

Global riverine loads in the context of published estimates
In our modelling framework, the global pre-industrial riverine exports to the ocean amount to 3.7 Tg P yr -1 , 27 Tg N yr -1 , 158 Tg Si yr -1 and 603 Tg C yr -1 , whereas 0.3 Tg P, 2.2 Tg N, 10 Tg Si and 19 Tg C are retained in endorheic catchments. The 30 fraction of Fe-P assumed to not desorb in the coastal ocean was also subtracted from the global loads (0.2 Pg P yr -1 ).
In the following, we compare the magnitudes of the modelled pre-industrial riverine exports of P, N, Si and C, as well as of their fractionations, with a wide range of estimates found in published literature (Table 3). This was done in comparison to 14  pre-industrial estimates directly, as well as to 1970 estimates by the NEWS2 study, in order to assess the modelled loads in the context of present-day estimations.  Meybeck and Vörösmarty (1999), 10 Resplandy et al. (2018), 11  , 12 Berner et al. (1983), 13 Amiotte Suchet and Probst (1995), 14 Mackenzie et al. (1998), 15 Cai (2011). The modelled global P loads are slightly below the P export estimate range of 4 -4.7 Tg P yr -1 reported in Compton et al. (2000), which was constructed by scaling measured catchment P loads of pristine rivers to the global river freshwater discharge (thus prehuman). A recent modelling study by Beusen et al. (2016), which takes P retention processes in rivers into account, suggests a lower global load of 2 Tg P yr -1 for the year 1900. The 1970 estimate (7.6 Tg P yr -1 ) provided by the NEWS2 study, which considers substantial anthropogenic P inputs, nevertheless suggests a steep 20th century increase in the global P 5 export to the ocean when considering all three pre-industrial estimates. The modelled DIP export to the ocean (0.5 Tg P yr -1 ) is situated at the top range of the prehuman estimates (0.3 -0.5 Tg yr -1 ) and well below estimates for the year 1970 (1.1 Tg P yr -1 ).

Species
A direct fractionation of the global P flux to DIP, DOP and POP is not provided in the Beusen et al. (2016) study. Furthermore, we estimate similar global loads of DOP to Compton et al. (2000) (around 0.1 and 0.2 Tg P yr -1 , respectively). The modelled DOP value is also much lower than the 1970 DOP estimate (0.6 Tg P yr -1 ), which is realistic due to the human-caused increase of DOP loads during the 20th century (Seitzinger et al., 2010). The modelled POP global load is larger than the estimate of Compton et al. (2000), which could be due to the POM C:P ratio of 122:1 chosen in our study. Strong uncertainties exist in the global C:P ratios for riverine POM, with Meybeck (1993) suggesting a mole ratio of around 56 C:P, whereas Ramirez and Rose (1992) estimate a ratio of around 499, which would strongly affect the results given here. The particulate P load given in 5 the NEWS2 study is vastly higher than the POP load modelled in our study, but a large fraction of the estimate is likely detrital particulate P, which we do not account for in the present study. The modelled Fe-P (1.0 Tg P yr -1 ) is slightly below the range estimated in Compton et al. (2000) (1.5 -3.0 Tg P yr -1 ). However, the assumed reactive fraction of the Fe-P loads (0.8 Tg P yr -1 ) here is close to how much Fe-P is suggested to be desorbed in the coastal ocean in Compton et al. (2000) (1.1-1.5 Tg P yr -1 ).

10
Despite our simplified assumption of N loads being directly coupled to P loads, the modelled global N load is situated within the prehuman and contemporary land-ocean N load estimates modelled in Green et al. (2004) (21 and 40 Tg N yr -1 , respectively). The modelled annual DIN load (3.4 Tg N) is slightly higher than the prehuman load given in the Green et al.
(2004) study (2.4 Tg N yr -1 ). In Beusen et al. (2016), the global pre-industrial N load is suggested to be lower (19 Tg N yr -1 ) due to in-stream retention and removal of N. 15 Since we assume that the change in the global DSi load over the 20th century is small, we directly compare our modelled pre-  Substantial amounts of particulate Si are suggested to be delivered to the ocean, yet it is not clear how much of it is dissolvable in the ocean (Tréguer and De La Rocha, 2013), and therefore, particulate Si is not considered in our study. The effects from contemporary increases in river damming on global DSi loads is also uncertain. Increased damming could have increased the global silica retention in rivers (Ittekkot et al., 2000;Maavara et al., 2014). Therefore, the global pre-industrial DSi load to the ocean could potentially have been larger than for the present day, but quantifying the implications of damming for the in-stream 25 retention of DSi, as well as of suggested increased DSi inputs from weathering and from the terrestrial biosphere (e.g. Gislason et al. (2009)) escape the scope of this study.
In our study, we also assumed that the pre-industrial C loads of all species were the same as for the present day. While the C retention along the land-ocean continuum might have increased during the 20 th century, enhanced soil erosion through changes in land use might have also increased the C inputs to the freshwater systems, leaving question marks on the magnitude of the 30 net anthropogenic perturbation Maavara et al., 2017). The modelled total C, DIC, DOC and POC fluxes agree with the present-day estimate ranges shown in Table 3, but are also situated at the lower end of these ranges. The organic C loads, which originate directly from the NEWS2 study, are particularly low.
The large spreads of the riverine export estimates underline difficulties in constraining pre-industrial riverine fluxes. Even for the present-day, Beusen et al. (2016) note large differences between the simulated loads of their study and the previous 35 global modelling study NEWS2. In turn, upscaling from point measurements are often based on data collected by Meybeck (1982) for the pre-1980 time period without taking into account more recent river measurement data. They also rely on the assumption of a linear relationship between riverine discharge and compound loads. While we acknowledge a certain degree of uncertainty in the numbers provided in this study, the modelling approach chosen nevertheless leads to global pre-industrial riverine loads that are in line with what was suggested previously and constitutes a framework that could be incorporated within 5 state-of-the-art Earth System Models.

Hotspots of riverine loads
We observe several regions of disproportionate contributions to global riverine loads (Figure 3, Table 4), with a distribution largely following spatial patterns of the weathering release. Our framework reveals large differences between the riverine exports from the Northern Hemisphere and from the Southern Hemisphere. The Northern Hemisphere thereby accounts for 10 an annual riverine C input of 404 Tg C to the ocean, vastly dominating the global loads (67% of global riverine C), while the Southern Hemisphere delivers 199 Tg C (33% global riverine C).
Rivers that drain into the Tropical Atlantic basin consist of a major fraction of the global oceanic biogeochemical supply (Table 4). This is due to major rivers of the South American continent debouching into the ocean basin, as well as considerable exports provided by West African rivers. According to our framework, the annual C loads of the seven largest rivers debouching and Niger, but are overestimated for the smaller catchments of the Paraíba, Volta and São Fransisco. Calculating the regional sum of the modelled pre-industrial and estimated present-day DIP loads reveals a much lower pre-industrial load (81.8 10 9 g P yr -1 ) with respect to the total present-day load (276 10 9 g P yr -1 ), suggesting a realistic increase of 194 10 9 g P yr -1 due to 25 anthropogenic perturbations.
Although less substantial in the context of global loads, major Arctic rivers (Yukon, Mackenzie, Ob, Lena, Yenisei) provide a large C supply to the Arctic Ocean, with a dominant fraction of the C load consisting of DIC. In our framework, these rivers provide 37.5 Tg DIC (10% of global DIC), 14.4 Tg DOC (11 % of global DOC) and 4.4 Tg POC annually to the Arctic Ocean.
The total riverine C loads to the Arctic therefore amount to 56 Tg C yr -1 , thus 9% of global C loads. The DIC, DOC and POC 5 Tg POC (Dittmar and Kattner, 2003). Total Arctic DIP loads (40.8 10 9 g P yr -1 ) derived from our modelling approach are slightly higher with regard to the published literature estimate of 35.8 10 9 g P yr -1 . DIP inputs from anthropogenic sources are considered to be small for Arctic catchments (Seitzinger et al., 2010), which explains why the modelled pre-industrial DIP loads are of similar magnitudes to observed DIP loads for the present day. Modelled Arctic DIC and DIP loads, which are strongly affected by inputs from the weathering models used in this study, are slightly larger with respect to published literature estimates, suggesting that the runoff scaling factor of 1.59 that we applied within the weathering models might be too high for this region. Certain factors, such as the effects of permafrost on weathering rates, which could impact the dynamics of weathering in the region, also remain enigmatic . 5 Southeast Asian rivers also provide large loads of biogeochemical tracers to the ocean. The Huang He, Brahmaputra-Ganges, Yangtze, Mekong, Irrawaddy and Salween rivers, which have catchment areas characterized by warm and humid climates, provide 92.4 Tg C yr -1 to the ocean (15% of global C loads). We observe vastly elevated DIP levels for the present-day estimates (76.4 10 9 g P yr -1 ) with regards to our modelled pre-industrial loads (299 10 9 g P yr -1 ). The NEWS2 study suggests a strong present-day perturbation of the DIP loads due to anthropogenic inputs to the region's catchments (Seitzinger et al., 4 Implications for the global ocean biogeochemistry

The ocean state -An increased biogeochemical coastal sink
In the ocean simulations, the magnitudes of oceanic nutrient inputs (P, N and Si) do not strongly differ between REF and RIV ( ocean therefore must act as an increased biogeochemical sink in RIV, since the produced organic matter reaches the shallow coastal sea floor faster than in the open ocean, which allows for less time for the organic matter to be remineralized within the Table 4. Regional hotspots of riverine C loads [Tg C yr -1 ] and DIP loads [10 9 g P yr -1 ] compared with regional estimates. 1 Araujo et al.  water column. This is confirmed by the increased organic matter flux to the global coastal sediment (defined as areas with less than 250m depths), with an increase from 0.18 Gt C yr -1 (REF) to 0.25 Gt C yr -1 (RIV). The global coastal POM deposition range given in a review by Krumins et al. (2013) shows a large degree of uncertainty (0.19-2.20 Gt C yr -1 ), but nevertheless hints that the coastal POM deposition flux is possibly improved in RIV.
While Alk inputs were also added at nearly the same levels in REF as in the calculated riverine loads of RIV, the global 5 C inputs are on the other hand increased by almost 100% in RIV in comparison to REF. The inorganic (366 Tg C yr -1 ) and organic (237 Tg C yr -1 ) C inputs in RIV show stronger agreement with the riverine inorganic (260-550 Tg C yr -1 ) and organic C (270 -350 Tg C yr -1 ) global load estimates found in literature (Meybeck, 1982;Amiotte Suchet and Probst, 1995;Mackenzie et al., 1998;Meybeck and Vörösmarty, 1999;Hartmann et al., 2009;Seitzinger et al., 2010;Cai, 2011;, signifying an improvement in the model's C inputs with respect to REF. These higher C inputs also result in a net long-term 10 outgassing flux (231 Tg C yr -1 ), which is suggested in literature for the pre-industrial time period and is absent in REF. The net outgassing largely originates from higher DIC:Alk ratios of the riverine inputs (1:1) than is exported through the net CaCO 3 production (1:2), as well as higher C:P ratios of tDOM inputs (2583:1) than is exported in the oceanic organic ratio (122:1).
Both of these imbalances between oceanic input and export ratios lead to an increase of the pCO 2 and CO 2 outgassing in the long-term (See Appendix C.2).
15  Boyer et al. (2013)), whereas the global mean surface DSi concentration is higher than in the global concentration derived from the WOA data set (13.6 and 7.5 µM DSi). The WOA dataset is constructed from present-day observations of an ocean state that might already be perturbed by a substantial increase in riverine P and especially N loads, whereas the model shows pre-industrial concentrations. Considering the stronger anthropogenic increase in DIN riverine loads (e.g. Seitzinger et al. (2010)) could plausibly shrink some of the disagreement with the WOA dataset in the case of oceanic DIN concentrations. A large part of the DIN underestimation in the model is however most likely due to notably large Tropical Pacific oxygen minimum zones, which cause a large DIN sink 5 due to denitrification. Furthermore, the lower surface concentrations of DIP and DIN than found in the WOA dataset suggest that the coastal sink of biogeochemical compounds might be too large, meaning that the exports of P and N to the open ocean might be too low in the model. Nevertheless, the surface DIN:DIP ratio in RIV is slightly improved in comparison to REF with respect to WOA data, which is most likely due to the shrinking of the Tropical Pacific oxygen minimum zone, induced by reduced nutrient supplies to the region, and the corresponding decrease in denitrification.

Riverine-induced NPP hotspots
The open ocean regions with the highest NPP rates in REF remain the most dominant areas of NPP in RIV (Figure 5a,b), which also confirms the major importance of the ocean physics in dictating spatial patterns of the global NPP. In comparison 20 to REF, substantial enhancements of the NPP are simulated near major river mouths in RIV. In proximity to river mouths in lower latitudes, the uptake of nutrients by the phytoplankton occurs efficiently due to favorable light conditions, thus adding to the local organic matter stock (Figure 5c,d). This can be observed on the shallow shelves of the western Tropical Atlantic, where major nutrient supplies are provided (see section 3.3). Moreover, the NPP is increased in certain semi-enclosed seas such as the Caribbean Sea, the Baltic Sea, the Black Sea and the Yellow Sea, where satellite observational data also suggest 25 high chlorophyll concentrations (Behrenfeld and Falkowski, 1997).  key reason explaining this inefficient export (Ichikawa and Beardsley, 2002). The decrease in the Equatorial Pacific NPP is also likely responsible for most of the shrinking of global oxygen minimum zone volume ( Table 5).
The Arctic Ocean does not undergo a noticeable increase in NPP, despite high nutrient concentrations in the basin in RIV.
This can be explained by the region's light limitation over most of the year, as well as the presence of sea ice coverage, which both inhibit the primary production especially during winter. In the entire basin, the nutrient concentrations are much higher 5 than what is suggested in the WOA database. In Bernard et al. (2011), where riverine nutrient inputs were added to the ocean according to the NEWS2 study, similarly high concentrations of DSi were found in the Arctic. Furthermore, Harrison and Cota (1991) suggest that in the Arctic Ocean, nutrients limit phytoplankton growth in the late Summer. Although the Arctic NPP in RIV and REF is substantially higher in the summer than in other seasons, the NPP is never nutrient-limited for the vast majority of the ocean basin.

Riverine-induced CO 2 Outgassing
The addition of riverine C loads causes a simulated net oceanic CO 2 source of 231 Tg C yr -1 to the atmosphere (Table 5). While the hotspots of the riverine-induced C outgassing are regions in proximity to major river mouths (Figure 6b  transfer of riverine C from the Northern Hemisphere to the Southern Hemisphere (Figure 6c). This interhemispheric transfer of C in the ocean has been a topic of discussion in literature, with studies of Aumont et al. (2001) and Resplandy et al. (2018) also suggesting a transport of C between latitudinal regions of the ocean that compensates for the heterogeneous distribution of riverine C inputs.
The Arctic rivers (Lena, Mackenzie, Yenisei, Ob, Oder, Yukon) also provide C to their respective shelves, which causes 10 a net outgassing of riverine C on the Laptev shelf and in the Beaufort Sea (2.2 Tg C yr -1 and 2.3 Tg C yr -1 , respectively).
The impacts of riverine C loads on the Arctic shelves can also be observed in the present-day coastal ocean pCO 2 data set of Laruelle et al. (2017), in which very high pCO 2 values are reported near river mouths (>400ppm).

Impacts of riverine loads on coastal region
With respect to the 10 shallow shelf regions chosen for this study (Table 1), the catchments of the low latitude regions (5. CSK, 6.BEN, 7.SEA, 8.TWA, 9.CG) provide substantially more C and nutrients to the coastal ocean than the high latitude regions in our framework (Figure 7), although the differing size of the coastal regions and of their catchments may play a strong role in explaining these differences. The tropical West Atlantic (8.TWA) has the largest riverine supplies of biogeochemical 5 compounds, which are largely provided by the Amazon and Orinoco rivers. In the tropical regions of the Bay of Bengal (6.BEN) and Southeast Asia (7.SEA), the fraction of C delivered as POC is substantially higher than for the rest of the regions (Figure 7b). Furthermore, in the high latitude regions (1.BS, 2.LS, 3.NS, 10.SAM), DIC loads are the major source of riverine C, whereas for the other regions, organic C is the largest contributor to the regional C load. The higher importance of DIC loads for the Arctic C supplies can also be observed in Table 4. 10 In RIV, major coastal ocean NPP increases of 166%, 377% and 71% for the tropical West Atlantic (3.TWA), Bay of Bengal (5.BEN) and East China Sea (6.CSK), respectively, are simulated, in comparison to REF (Figure 8b). The availability of light, as well as the large riverine supplies of nutrients to these regions (Table 4 in Section 3.2), provide optimal conditions to enhance the biological production. Surprisingly, the Southeast Asian shelf (6.SEA) does not undergo a large relative increase (4.OKH) regions substantial NPP increases due to the addition of riverine inputs are also not simulated. Although the NPP is strongly enhanced in the direct proximity to the Paraná river (Figure 5b), the vastness of the South American shelf (1553 10 9 m 2 ) also makes the region less sensitive to riverine inputs. In published literature, the nutrient supply which drives the NPP on the Patagonian Shelf is also confirmed to be strongly controlled by open ocean inflows, and river supplies are suggested to only have a limited influence in the region (Song et al., 2016). 15 The Arctic shelf regions also do not undergo a strong NPP increase in RIV with respect to REF (8% and 5% increases for the Beaufort Sea, 1.BS, and Laptev Sea, 2.LS, respectively). We however do not consider a seasonality of the riverine inputs.
Larger inputs of nutrients in months of larger discharge from April to June (Le Fouest et al., 2013) could cause a more efficient usage of the riverine nutrients on the shelves, since the light availability is higher and the sea ice coverage is reduced during these months. 20

27
An increase in CO 2 outgassing due to the increased C inputs is simulated in all regions in RIV. In the Arctic regions (Beaufort Sea and Laptev Sea), the relative change in the air-sea CO 2 flux is very pronounced, whereas the impact is generally not as strong in the lower latitude regions due to the enhancement of the biological DIC uptake caused by riverine nutrient inputs.
The tropical West Atlantic is an exception to this latitudinal pattern, since the large riverine C supplies also cause a substantial change in the CO 2 flux of the region. In the North Sea, we observe an enhancement of C outgassing induced by riverine loads, 5 but the region remains a substantial pre-industrial sink of atmospheric CO 2 in RIV, as is still suggested for the present day by 6 Origins and fate of riverine carbon In our simplified model framework (Figure 9), we quantify the land sources of riverine C (1-3), its riverine transfer to the ocean (4), and the long-term fate of the riverine C in the ocean (5)(6)(7)(8)(9)(10). While the terrestrial fluxes are derived from the weathering and organic matter export models, the fluxes of the riverine C in the ocean are derived from RIV-REF differences of the ocean simulation. The long-term oceanic CO 2 flux is furthermore decomposed to illustrate the contributions of riverine inorganic and 5 organic C inputs to the oceanic outgassing flux in a model equilibrium analysis. The detailed calculations of the land and ocean fluxes are shown and explained in detail in the Appendix C (C1 for terrestrial and C2 for oceanic fluxes).
In our framework, the global pre-industrial terrestrial uptake of atmospheric C and its export to rivers as DIC amounts to 529 Tg C yr -1 . The sink consists of 280 Tg C yr -1 from the drawdown induced by weathering (1) and of 249 Tg C yr -1 from the terrestrial biological uptake (3). During the weathering process, 94 Tg C yr -1 is moreover released from the lithology through 10 carbonate weathering (2). During silicate weathering, all the C released originates from atmospheric CO 2 . The terrestrial biological uptake (3) is derived from the net global export of organic C to the ocean. It therefore implicitly takes into account the net land biological uptake and its export to freshwaters, as well as all net sinks and sources in river systems. A total of 603 Tg C yr -1 is transferred laterally to the ocean (4) while taking into consideration an endorheic catchment loss of 19 Tg C yr -1 .
In the ocean, riverine inputs of C cause a long-term net annual atmospheric source of 231 Tg C yr -1 (Table 5). We propose 15 a decomposition of this long-term CO 2 flux into sources and sinks induced by the inputs of the riverine compounds (Appendix C.2). Assuming model equilibrium, the oceanic outgassing flux can be decomposed into a source from riverine inorganic C inputs (183 Tg C yr -1 , 5), a source from terrestrial organic C (128 Tg C yr -1 ), 6) and a sink caused by the enhancement of the biological production prompted by the addition of dissolved inorganic nutrients and its associated alkalinity production (69 Tg C yr -1 , 7). We also calculate a sink due to disequilibrium at the atmosphere-water column interface in the model (11 20 Tg C yr -1 , Dr1). The production and the downwards export of CaCO 3 and POM within the ocean lead to simulated sediment deposition fluxes of 188 Tg C yr -1 for CaCO 3 (inorg. C. flux, 8) and of 582 Tg C yr -1 for POM (org. C flux, 9). The dissolution of CaCO 3 and the remineralization of POM within the sediment lead to a DIC flux of 385 Tg C yr -1 from the sediment pore water to the water column. The net C flux at the sediment interface (8+9-10) is therefore a burial flux of 385 Tg C yr -1 . The calculated equilibrium C burial flux, which is the difference between the riverine C inputs and the equilibrium CO 2 outgassing 25 (4-5-6+7), is 361 Tg C yr -1 , which implies that there is a deviation of 24 Tg C yr -1 (Dr2) towards the sediment. with the simulated C burial (385 Tg C yr -1 ) deviating from the calculated equilibrium state C burial (361 Tg C yr -1 ). The similar deviations from the equilibrium states at the atmosphere-water column and water column-sediment interfaces suggest that the model disequilibrium, which is likely due to long time scales of the processes taking place in the sediment, translates relatively efficiently into a perturbation of the CO 2 flux at the atmosphere-water column interface. 30 The simulated riverine-induced oceanic CO 2 outgassing flux of 231 Tg C yr -1 is consistent with the estimate range of 200-400 Tg C yr -1 given in Sarmiento and Sundquist (1992), who derived an annual global riverine C load of 300-500 Tg C, and with Jacobson et al. (2007) and Gruber et al. (2009), who suggest a slightly higher natural CO 2 outgassing flux of 450 Tg C Net land biological C uptake (derived directly from riverine organic C exports). 4. Total riverine C exports. 5. Oceanic C outgassing from riverine DIC (i.e. inorganic C) inputs. 6. Oceanic C outgassing resulting from riverine tDOM inputs.
7. Oceanic C uptake due to the enhanced primary production by dissolved inorganic nutrients and the corresponding alkalinity production. 8.  Resplandy et al. (2018) proposes riverine C loads and oceanic C outgassing of both 780 Tg C yr -1 . It is however unclear if and how the oceanic carbon removal of the riverine C through sediment burial was considered in the study.
Furthermore, we observe an imbalance in the calculated pre-industrial land uptake and the oceanic outgassing of atmospheric CO 2 in our approach, with the land uptake outweighing the oceanic outgassing (combined net atmospheric sink of 299 Tg C yr -1 ). Accounting for further sources of atmospheric CO 2 such as volcanic emissions and shale oxidation would therefore be 5 necessary to achieve a stable atmospheric carbon budget in a fully coupled land-ocean-atmosphere setting, since ESMs assume constant pre-industrial atmospheric CO 2 levels. For instance, Mörner and Etiope (2002) suggest long-term volcanic annual emissions in the range of 80-160 Tg C yr -1 , whereas Burton et al. (2013) estimate a long-term volcanic outgassing flux as high as that of silicate weathering drawdown. In our approach, the modelled silicate weathering CO 2 drawdown is 196 Tg C yr -1 .
Additionally to the volcanic CO 2 emissions, the global atmospheric CO 2 land source from the oxidation of organic carbon 10 in shale of around 100 Tg C yr -1 given by Sarmiento and Sundquist (1992) would then approximately close the atmospheric carbon budget in our framework.

Rivers in an Earth System Model setting
Our approach, which quantifies riverine exports as a function of the climate variables precipitation, surface runoff and temperature, can potentially be used to estimate land-ocean exports in an ESM setting. Since we establish a baseline for the pre-industrial state of the ocean while considering riverine loads, our study could also be used to assess oceanic impacts of 5 perturbations of riverine loads due to temporal changes in weathering rates (Gislason et al., 2009), or due to increased anthropogenic P and N inputs to catchments (Seitzinger et al., 2010;Beusen et al., 2016).
This study however also relies on strong assumptions in order to perform simulations at the global scale, although improvements within the framework are possible. For one, the weathering mechanisms, the processing uptake of nutrients by the terrestrial biology, the representation of hydrological flow characteristics, as well as transformation processes and retention in

Dynamics of terrestrial organic matter in the ocean
The composition of terrestrial organic matter, as well as its fate in the ocean have been strongly debated in the past. Recent 20 work shows that despite its low biological reactivity, tDOM can be mineralized by abiotic processes such as photodegradation in the ocean (Fichot and Benner, 2014;Müller et al., 2016;Aarnos et al., 2018), despite already having been strongly degraded during its transit in rivers. The global magnitude of the degradation of tDOM in the ocean is however strongly uncertain. Few studies in the past have addressed the composition of POM, although it is also thought to be efficiently remineralized in the coastal ocean sediment (Hedges et al., 1997;Cai, 2011). While the C loads from POM are the lowest loads of all C compounds 25 considered in this study, a differing C:P ratio to the one chosen in this study could affect the modelled riverine-induced C outgassing flux in the ocean.
Rates of remineralization processes have been suggested to be higher in the coastal ocean than in the open ocean (e.g. Krumins et al. (2013)). Using a higher coastal organic matter remineralization rate in the sediment module of HAMOCC could potentially reduce the coastal biogeochemical sink described in this study and increase the offshore export of biogeochemical 30 compounds.

Arctic Ocean
Simulated nutrient concentrations are particularly high in the Arctic Ocean with regard to WOA data, suggesting that this region with strong riverine inputs might be poorly represented in the HAMOCC. Difficulties to represent the region could be due to fine circulation features. For instance, outflows through narrow passages have been shown to be affected by model resolution (Aksenov et al., 2010). The simulated sea ice coverage in our study, which is around 85% for the Laptev Sea during the whole 5 year on average, could be overestimated. Moreover, the NPP in the region might be underestimated due to photosynthesis taking place under ice, in ice ponds and over extended daytime periods in the summer months (Deal et al., 2011;Sørensen et al., 2017), all of which are not represented in HAMOCC.

Summary and conclusions
In this study, we account for weathering and non-weathering inputs to river catchments to quantify global annual pre-industrial 10 loads of 3.7 Tg P, 27 N, 168 Tg DSi, and 603 Tg C to the ocean. These loads are consistent with published literature estimates, although we acknowledge a certain degree of uncertainty regarding their magnitudes. While we omit the in-stream retention of P during its riverine transport, which reduces the global P exports to the ocean (Beusen et al., 2016), our estimate of the global pre-industrial P export to the ocean is comparable in magnitude to an estimation that determines riverine P exports by upscaling from pristine river measurements (Compton et al., 2000). 15 We identify the Tropical Atlantic catchments, the Arctic Ocean, Southeast Asia and Indo-Pacific islands as regions of dominant contributions of riverine supplies to the ocean. These 4 regions account for over 51% of land-ocean C exports in total, with Tropical Atlantic catchments supplying around 20% of C to the ocean globally. We also observe that the contributions of different C species differ between the regions. Most prominently, the C supply of the Indo-Pacific islands is dominated by particulate organic C loads, which have been identified to be more strongly controlled by extreme hydrological events than 20 other C species (Hilton et al., 2008).
In the ocean model, riverine inputs of C lead to a net global oceanic outgassing flux of 231 Gt C yr -1 , a comparable value with regards to previous estimates of 200-450 Tg C yr -1 (Sarmiento and Sundquist, 1992;Jacobson et al., 2007;Gruber et al., 2009).
This outgassing flux can be decomposed into two source terms caused by inorganic C inputs (183 Tg C yr -1 ) and organic C inputs (128 Tg C yr -1 ), and a net sink term (80 Tg C yr -1 ) caused by the enhanced biological C uptake due to riverine inorganic 25 nutrient supplies, corresponding alkalinity production and a slight model drift in alkalinity. The magnitude of the outgassing is however strongly dependent on the magnitude of riverine C loads. In our framework, both riverine C loads and their induced oceanic C outgassing are thereby situated within the lower range of estimates.
We find evidence of a substantial oceanic interhemispheric transport of riverine C from the Northern Hemisphere to the Southern Hemisphere, with a larger relative C outgassing flux simulated over the Southern Hemisphere (49% of global out-30 gassing) than its relative riverine C loads (36% of global C loads). We also show that the Southern Ocean outgasses 17 Tg of riverine C, despite not having a direct riverine source of C, meaning that riverine C is transported within the ocean interior to the Southern Ocean in the model. This interhemispheric transfer of riverine C in the ocean has previously been suggested to contribute to the pre-industrial Southern Ocean source of atmospheric CO 2 for the pre-industrial time frame (Sarmiento et al., 2000;Aumont et al., 2001;Gruber et al., 2009;Resplandy et al., 2018). Here, we show that riverine C exports derived from state-of-the-art land export models confirm the larger contribution of the Northern Hemispheric terrestrial C supply to the ocean. Part of this uneven terrestrial C supply is then compensated by the transport of C within the ocean and is outgassed remotely.

5
Our results help identify oceanic regions that are sensitive to riverine fluxes. Riverine-induced changes in the regional NPP are mostly found in coastal regions, but substantial riverine-caused CO 2 outgassing can also be observed in the open ocean of the Tropical Atlantic. In general, we show latitudinal differences in the sensitivity of the NPP and of the CO 2 fluxes to the riverine inputs in various coastal regions. While a high sensitivity in the NPP is found in the tropics and subtropics, with the tropical West Atlantic, the Bay of Bengal and the East China Sea showing large NPP increases of 166%, 377% and 71%, 10 respectively, the relative changes in the regional CO 2 fluxes are mostly larger at higher latitudes. For instance, the Laptev Sea and the Bay of Beaufort become atmospheric sources of 2.2 Tg C yr -1 and 2.3 Tg C yr -1 respectively, despite being sinks of atmospheric CO 2 without accounting for riverine inputs. While our analysis revolves around the implications of preindustrial riverine loads, regions that show strong responses to the riverine loads in this study might also be strongly affected by anthropogenic perturbations of riverine exports. 15 Deriving riverine exports as a function of ESM variables enables a representation of the riverine loop; starting from the terrestrial uptake of C, to its riverine export and ending with its long-term outgassing in the ocean or export to the oceanic sediment. In order to implement the framework into a coupled land-atmosphere-ocean setting such as an ESM, the atmospheric pre-industrial budget would have to be balanced. We emphasize the need to consider long-term terrestrial CO 2 sources originating from long-term volcanic activity and from shale oxidation in order to compensate for imbalances between the terrestrial 20 C sink and the oceanic C source in our framework.
Throughout this study, we find global heterogeneity in the spatial features of weathering fluxes, riverine loads and their implications for the ocean biogeochemistry. We confirm that considering riverine exports to the ocean is central when assessing the biogeochemistry of coastal regions, but also find implications for certain open ocean regions (i.e. the Tropical Atlantic).
Our study also shows the necessity to account for the river-induced oceanic C outgassing in ocean biogeochemistry models, 25 since our conservative estimate consists of around 10% of the magnitude of the present-day oceanic C uptake.
Code and data availability. Code, primary data and scripts needed to reproduce the analyses presented in this study are archived by the Max Appendix B: Runoff, temperature and precipitation The weathering yields provided in this study are dependent on the MPI-ESM spatial representation of pre-industrial surface runoff, surface temperature and precipitation ( Figure B1a,b,c). Previous work by Goll et al. (2014) describes the MPI-ESM performing well when estimating surface temperatures at single 5 grid cells with regard to observations. The global precipitation in MPI-ESM is slightly higher than the global precipitation that is reported in the Precipitation Climatology Project (GCPC) (Adler et al., 2003), which is discussed along with spatial biases of the precipitation in Stevens et al. (2013). Most notably, the precipitation is too strong over extratropical land surface and too little over tropical land surface. Runoff is less well reproduced globally. For the given time period of the simulation used in this study, the global runoff is 23,496 km 3 yr -1 . This is significantly lower than the global runoff estimations of 36,600-38,300 km 3 yr -1 (Fekete et al., 2002;Dai and Trenberth, 2002). The difficulty of representing several processes that control the runoff, 5 such as evapotranspiration and condensation, is also reflected in the large range of global runoff means of other Earth System Models, which range from 23,000 to 42,500 km 3 yr -1 (Goll et al., 2014). The spatial patterns in the CMIP5 simulation are however comparable with the annual runoff patterns reported in Fekete et al. (2002), with high surface runoff observed in the Amazon basin, West Africa, Indo-Pacific Islands, Southeast Asia, eastern North America, Northern Europe as well as in Siberia. Due to the strong underestimation of the model regarding the runoff in relation to the combined runoff mean of Fekete 10 et al. (2002) and Dai and Trenberth (2002), we conclude that a scaling factor of 1.59 is necessary to represent runoff plausibly at the global scale. The global runoff from OMIP, which provides freshwater to the ocean model, is however more plausible (32,542 km 3 yr -1 ) and was therefore not scaled. We scaled the model runoff to the global runoff estimation from Fekete et al. (2002) and not to the OMIP runoff, since scaling to the OMIP runoff still led to lower global P and DIC/Alk release values than were reported in literature estimates.

15
Appendix C: Derivation of carbon fluxes in the simplified coupled system C1 Terrestrial fluxes 1. Carbonate and silicate weathering cause a terrestrial uptake of 280 Tg C yr -1 from the atmosphere according to the weathering model simulations (Table 2).

2.
Carbonate mineral weathering causes a lithological C release flux of 94 Tg C yr -1 DIC, as shown in Section 3.2. 20 3. Carbon contained in terrestrial organic matter exports originates from the atmosphere (Meybeck and Vörösmarty, 1999).
The net C uptake by the terrestrial and riverine biology is therefore the same as the lateral organic C export, which we derived from NEWS2 (Seitzinger et al., 2010) DOC and POC exports. The net atmospheric C uptake by the terrestrial biology, while taking into account all respiration processes on land and in rivers, is the sum of the lateral POC and DOC exports (249 Tg C yr -1 ).

4.
The global riverine C export to the ocean consists of the sum of weathering and organic matter C exports (623 Tg C yr -1 ), minus a loss term due to endorheic river retention (19 Tg C yr -1 ). This results in the global riverine C export of 603 Tg C yr -1 (values are rounded).

C2 Long-term ocean fluxes
In the model, riverine loads cause oceanic C outgassing due to the inputs of inorganic C and tDOM, while the inputs of dissolved inorganic nutrients cause a sink of atmospheric C through their enhancement of the biological C uptake. Increased primary production thereby also produces alkalinity, which also potentially causes an uptake of atmospheric C.
Inorganic C is delivered to the ocean by rivers as 1 mol DIC and 1 mol alkalinity, which represents the supply of HCO 3 -, 5 and exported from the ocean as CaCO 3 , which accounts for 0.5 mol DIC per 1 mol alkalinity, leaving a surplus of 0.5 mol DIC and 0 mol alkalinity when balancing over the oceanic inputs and exports. Increasing the oceanic DIC pool without increasing the alkalinity directly increases the dissolved CO 2 concentrations, which in turn causes outgassing: Equilibrium model outgassing caused by riverine inorganic C inputs is therefore 0.5 fold of the riverine DIC loads. 10 Outgassing from organic matter inputs results from the high C:P ratio of tDOM. tDOM, in contrary to POM, is not exported to the sediment. It is mineralized in the ocean providing dissolved inorganic compounds in the C:N:P ratio of 2583:16:1, but the subsequent uptake of the released inorganic compounds happens at a C:N:P ratio of 122:16:1, resulting in a DIC overshoot.
Since the net alkalinity over the entirety of Eq. C2 is also constant, the DIC increase causes a pCO 2 increase (*simplified equation):
The riverine loads of DIP and DIN on the other hand cause C uptake through their enhancement of the biological primary production. DIC/CO 2 is thereby removed, thus sinking pCO 2 (*simplified equation): DIP + 16DIN + * +122CO 2 => C 122 P N *

16
(C3) 20 The resulting C uptake from the equation is therefore 122-fold the mole DIP inputs. Additionally, alkalinity is produced in Eq. C3. The uptake of DIN (assuming nitrate and not ammonium,as is done in HAMOCC) and DIP through primary production causes a net alkalinity increase by a factor of Alk:P = 17:1 (Wolf-Gladrow et al., 2007). In model equilibrium, the produced alkalinity must also be exported through Eq. C1. This alkalinity is exported in a C:Alk ratio of 1:2, meaning an uptake of 0.5 DIC per mol alkalinity exported. The C uptake enhancement from the alkalinity increase is therefore the 17 * 0.5 -fold of the 25 (bioavailable) DIP loads.
With the help of these equations, the fluxes of riverine C in the ocean of Figure 9 were calculated: 37 5. According to the CaCO 3 export equation (Eq. C1), half of the riverine DIC input is exported to the sediment as CaCO 3 and the other half is outgassed as CO 2 in model equilibrium state. The contribution of outgassing caused by inorganic C in the ocean is therefore 0.5-fold the DIC inputs (366 Tg C yr -1 ) and therefore 183 Tg C yr -1 assuming model equilibrium.
6. tDOM input C:P rations vastly exceed the oceanic sediment export C:P ratios of organic matter, which causes model equilibrium outgassing in the ocean. Eq. C2 shows that for every mol tDOM supplied to the ocean, in model equilibrium 5 122/2583 of C is exported to the sediment and 2461/2583 of C increases the dissolved CO 2 pool, which is outgassed in the long term. The equilibrium outgassing is therefore the tDOM carbon load (134 Tg C yr -1 ) multiplied by 2461/2583, which results in 128 Tg C yr -1 . In the case of POM, since the C:P ratio of the riverine input (122:1) is the same as the ratio of the export to the sediment in the ocean, there is no effect on the long-term equilibrium C outgassing flux.

7.
Since P has no further sinks or sources in the model other than riverine inputs and sediment burial as organic matter, in 10 model equilibrium, the same amount of P supplied by rivers is buried in the sediment. When DIP is taken up by the biology and transformed to organic matter, carbon is also taken up in a mole ratio of C:P = 122:1. Accounting for this uptake through the biological production enhancement by DIP inputs (including all bioavailable DIP) from rivers (1.4 Tg P yr -1 ) results in the uptake of 67 Tg C yr -1 through the biological enhancement. Furthermore, when DIP and DIN are transformed to organic matter as organic matter, an alkalinity increase of 17 mol per mol DIP uptake takes place (Eq. C3, Wolf-Gladrow et al. (2007)). This 15 increase in alkalinity causes the further uptake uptake and export of 2 Tg C yr -1 , resulting in a total sink of 69 Tg C yr -1 .
Dr1. We attribute the difference between the equilibrium CO 2 flux of 242 Tg C yr -1 and the modelled net CO 2 flux of 231 Tg C yr -1 (=11 Tg C yr -1 ) to the small surface alkalinity increase in the model simulations due to slight disequilibrium over the analysis time period.
8. The simulated global particulate inorganic C sediment deposition flux in HAMOCC is 188 Tg C yr -1 . 20 9. The simulated global organic C sediment deposition flux in HAMOCC model is 582 Tg C yr -1 .

10.
The modelled global DIC flux from the sediment back to the water, which originates from DIC release due to POM remineralization and CaCO 3 dissolution in the sediment, is 385 Tg C yr -1 .
Dr2. In model equilibrium, the net sediment burial C flux is the total riverine C inputs of 603 Tg C yr -1 (4) subtracted by the equilibrium outgassing of 242 Tg C yr -1 (5+6-7), which results in 361 Tg C yr -1 . The simulated net C burial flux in the 25 simulations of 385 Tg C yr -1 (8+9-10) deviates from the calculated model burial equilibrium flux. Therefore, the drift at the sediment-ocean interface is 24 Tg C yr -1 (385 Tg C yr -1 -361 Tg C yr -1 ), and is in the direction of the sediment.   Gaillardet, J., Dupré, B., Louvat, P., and Allègre, C.: Global silicate weathering and CO2 consumption rates deduced from the chemistry of