2.1 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 directions of the model at a horizontal resolution of 0.5^∘. 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.

We considered both weathering as well as non-weathering sources of P, N, Si, Fe, C and Alk to river catchments (Fig. 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₂ 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 ocean in our approach (Fig. 1).

The surface runoff, surface temperature and precipitation data used to drive the weathering release and land export models were obtained from the 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^∘. The runoff data were 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 that is discussed in Appendix B.

2.1.1 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 modeled anthropogenic perturbations to riverine loads for 1970 and 2000, only simulated small changes in global POC and DOC loads over this time period. Regnier et al. (2013) suggest a total anthropogenic perturbation 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 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 the ocean biogeochemistry model, which was also the reasoning for choosing the P:Fe mole ratio of $\mathrm{1}:\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{4}}$. The $\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{Fe}$ mole ratio of tDOM was therefore $\mathrm{2583}:\mathrm{16}:\mathrm{1}:\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{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 $\mathrm{C}:\mathrm{P}=\mathrm{122}:\mathrm{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 $\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{Fe}$ ratio as of oceanic POM ($\mathrm{122}:\mathrm{16}:\mathrm{1}:\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{4}}$).

Figure 1(a) Table of sources of nutrient, carbon and alkalinity inputs to river catchments and (b) scheme of origins and transformations of the catchment compounds, as well as biogeochemical processes in the ocean in our framework. The abbreviations are: Inorg. Comp.: Inorganic Compounds, tDOM: terrestrial dissolved organic matter, POM: particulate organic matter, DIP: dissolved inorganic phosphorus, Fe-P: iron-bound phosphorus, DIN: dissolved inorganic nitrogen, DSi: dissolved silica, DFe: dissolved iron, DIC: dissolved inorganic carbon, Alk: alkalinity.

Download

2.2 Ocean model setup

2.3 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 riverine loads that cover a variety of latitudes (Table 1). The shelves were defined to have depths of less than 250 m. The cutoff sections perpendicular to the coast were done according to MARgins and CATchments 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.

Table 1Comparison of the surface areas (10⁹ m²) of selected coastal regions with depths of under 250 m in the ocean model setup and in the segmentation approaches. The comparisons were done with the MARCATS (Laruelle et al., 2013) or COSCATs (Beaufort Sea, Laptev Sea, North Sea, Congo shelf, Meybeck et al., 2006). The shelf classes were defined as in Laruelle et al. (2013).

Download Print Version | Download XLSX

3.1 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 their spatial distributions (Fig. 2). The calculated global P release is 1.34 Tg P yr⁻¹ when considering the runoff scaling factor of 1.59, which fits within the P release range of 0.8–4 Tg P yr⁻¹ 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 Goll et al. (2014) estimates (1.1 and 0.8–1.2 Tg P yr⁻¹, respectively) originate from the same weathering model as is used here, the 1.9 Tg P yr⁻¹ 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 land–ocean 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 runoff estimates, thereby yielding an estimate of 2.5–4.0 Tg P yr⁻¹. Both estimates that originate from the upscaling of river measurements are therefore higher than the P weathering fluxes provided in the modeling 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⁻¹. This is within the range estimated by Beusen et al. (2009) (158–199 Tg Si yr⁻¹), who used the same model while using present-day observational data to drive the model for precipitation. Our estimate is also comparable to the 173 Tg Si yr⁻¹ value reported by Dürr et al. (2011).

The modeled DIC and Alk release amount to 374 Tg C yr⁻¹ and 37.8 Teq yr⁻¹ globally, which also take 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⁻¹ and originates directly from weathering. Further modeling studies also provide similar weathering release estimates of 260 to 300 Tg C yr⁻¹ (Berner et al., 1983; Amiotte Suchet and Probst, 1995). The atmospheric CO₂ drawdown induced by weathering is directly related to the weathering release of DIC, since silicate weathering draws down 1 mole of CO₂ per mole ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and carbonate weathering draws down 0.5 mol of CO₂ per mole ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ released (Eqs. 9 and 10). While we model a CO₂ drawdown flux of 280 Tg C yr⁻¹ induced by weathering, drawdown fluxes of 220–440 Tg C yr⁻¹ are reported in published literature (Gaillardet et al., 1999; Amiotte Suchet et al., 2003; Hartmann et al., 2009; Goll et al., 2014). The results imply that, of the 374 Tg C yr⁻¹ DIC released by weathering, 280 Tg C yr⁻¹ originates from atmospheric C drawdown, while the rest originates from the weathering DIC release from the carbonate lithology (94 Tg C yr⁻¹).

Table 2Global pre-industrial weathering release of P, Si, DIC and Alk, as well as CO₂ drawdown, modeled in this study and previously estimated in published literature.

Download Print Version | Download XLSX

Figure 2Weathering release yields of (a) P (kg P km⁻²), (b) Si (t Si km⁻²) and (c) DIC (t C km⁻²).

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 (Fig. 2). The disproportionate importance 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).

3.2 Global riverine loads in the context of published estimates

In our modeling framework, the global pre-industrial riverine exports to the ocean amount to 3.7 Tg P yr⁻¹, 27 Tg N yr⁻¹, 158 Tg Si yr⁻¹ and 603 Tg C yr⁻¹, whereas 0.3 Tg P, 2.2 Tg N, 10 Tg Si and 19 Tg C are retained in endorheic catchments. The fraction of Fe-P assumed to not desorb in the coastal ocean was also subtracted from the global loads (0.2 Pg P yr⁻¹).

Table 3Comparison of modeled global riverine loads (Model. global loads) with previous estimates (Tg yr⁻¹). The total loads thereby exclude particulate inorganic loads, except for the 1970 POP estimate, which includes all particulate P compounds from the NEWS2 study. Fe-P only considers all particulate P (PP) that is assumed to be desorbed in the ocean. The PP value from the NEWS2 study considers both POP as well as particulate inorganic phosphorus.

¹ Compton et al. (2000). ² Beusen et al. (2016). ³ NEWS2 (Seitzinger et al., 2010). ⁴ Beusen et al. (2009). ⁵ Dürr et al. (2011). ⁶ Tréguer and De La Rocha (2013). ⁷ Green et al. (2004). ⁸ Jacobson et al. (2007). ⁹ Meybeck and Vörösmarty (1999). ¹⁰ Resplandy et al. (2018). ¹¹ Regnier et al. (2013). ¹² Berner et al. (1983). ¹³ Amiotte Suchet and Probst (1995). ¹⁴ Mackenzie et al. (1998). ¹⁵ Cai (2011).

Download Print Version | Download XLSX

In the following, we compare the magnitudes of the modeled 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 pre-industrial estimates directly, as well as to 1970 estimates by the NEWS2 study, in order to assess the modeled loads in the context of present-day estimations.

The modeled global P loads are slightly below the P export estimate range of 4–4.7 Tg P yr⁻¹ 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 modeling 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⁻¹ for the year 1900. The 1970 estimate (7.6 Tg P yr⁻¹) provided by the NEWS2 study, which considers substantial anthropogenic P inputs, nevertheless suggests a steep 20th century increase in the global P export to the ocean when considering all three pre-industrial estimates. The modeled DIP export to the ocean (0.5 Tg P yr⁻¹) is situated at the top range of the prehuman estimates (0.3–0.5 Tg  yr⁻¹) and well below estimates for the year 1970 (1.1 Tg P yr⁻¹). 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⁻¹, respectively). The modeled DOP value is also much lower than the 1970 DOP estimate (0.6 Tg P yr⁻¹), which is realistic due to the human-caused increase of DOP loads during the 20th century (Seitzinger et al., 2010). The modeled 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 the NEWS2 study is vastly higher than the POP load modeled 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 modeled Fe-P (1.0 Tg P yr⁻¹) is slightly below the range estimated in Compton et al. (2000) (1.5–3.0 Tg P yr⁻¹). However, the assumed reactive fraction of the Fe-P loads (0.8 Tg P yr⁻¹) 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⁻¹).

Despite our simplified assumption of N loads being directly coupled to P loads, the modeled global N load is situated within the prehuman and contemporary land–ocean N load estimates modeled in Green et al. (2004) (21 and 40 Tg N yr⁻¹, respectively). The modeled 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⁻¹). In Beusen et al. (2016), the global pre-industrial N load is suggested to be lower (19 Tg N yr⁻¹) due to in-stream retention and removal of N.

Since we assume that the change in the global DSi load over the 20th century is small, we directly compare our modeled pre-industrial DSi load estimate with present-day DSi load estimates from published literature. The modeled global load of DSi is 158 Tg Si yr⁻¹ in our framework, which is situated at the low end of the present-day estimate range of 158–200 Tg Si yr⁻¹. The NEWS2 study used the same Beusen et al. (2009) DSi export model forced with present-day observational precipitation data. Dürr et al. (2011) and Tréguer and De La Rocha (2013) upscaled from discharge-weighted DSi concentrations at river mouths. 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 are 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 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 20th 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 net anthropogenic perturbation (Regnier et al., 2013; Maavara et al., 2017). The modeled 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 global modeling 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 modeling 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 state-of-the-art Earth system models.

3.3 Hotspots of riverine loads

We observe several regions of disproportionate contributions to global riverine loads (Fig. 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 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 in the basin (Orinoco, Amazon, São Francisco, Paraíba do Sul, Volta, Niger and Congo) amount to 123 Tg C (fractionated into 58 Tg DIC, 44 Tg DOC, 21 Tg POC), which consists of around 20 % with respect to global riverine C exports. These regional C loads agree very well with estimated values derived from monthly data of Araujo et al. (2014) (53 Tg DIC, 46 Tg DOC). Of these major rivers, the pre-industrial Amazon river provides the largest input of biogeochemical tracers to the ocean (modeled annual loads of 0.07 Tg DIP, 0.5 Tg DIN, 15.2 Tg DSi, 33.2 Tg DIC, 28.3 Tg DOC Tg, 17.1 Tg POC), which compare well with the present-day data (0.22 Tg DIP, 17.8 Tg DSi, 32.7 Tg DIC, 29.1 Tg DOC, Araujo et al., 2014). For the other catchments debouching in the tropical Atlantic, the modeled DIC loads are close to estimated values for the Congo, Orinoco and Niger, but are overestimated for the smaller catchments of the Paraíba, Volta and São Francisco. Calculating the regional sum of the modeled pre-industrial and estimated present-day DIP loads reveals a much lower pre-industrial load (81.8×10⁹ g P yr⁻¹) with respect to the total present-day load (276×10⁹ g P yr⁻¹), suggesting a realistic increase of 194×10⁹ g P yr⁻¹ due to 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⁻¹, thus 9 % of global C loads. The DIC, DOC and POC load levels are comparable to estimates of 29 Tg DIC yr⁻¹ (Tank et al., 2012), 17 Tg DOC yr⁻¹ (Raymond et al., 2007) and 5 Tg POC (Dittmar and Kattner, 2003). Total Arctic DIP loads (40.8×10⁹ g P yr⁻¹) derived from our modeling approach are slightly higher with regard to the published literature estimate of 35.8×10⁹ g P yr⁻¹. DIP inputs from anthropogenic sources are considered to be small for Arctic catchments (Seitzinger et al., 2010), which explains why the modeled pre-industrial DIP loads are of similar magnitudes to observed DIP loads for the present day. Modeled 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 (Hartmann et al., 2014).

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⁻¹ to the ocean (15 % of global C loads). We observe vastly elevated DIP levels for the present-day estimates (76.4×10⁹ g P yr⁻¹) with regards to our modeled pre-industrial loads (299×10⁹ g P yr⁻¹). 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., 2010), which can plausibly explain these differences.

The Indo-Pacific islands have been identified as a region with much higher weathering yields than the global average (Hartmann et al., 2014). Although this region only accounts for around 2 % of the global land surface, it provides 7 % (39 Tg) of C and most notably 10 % of the global POC (10 Tg C) delivered by rivers to the ocean annually, making the region a stronger terrestrial source of POC than all Arctic catchments. Thus, the POM mobilization through soil erosion is likely a major driver of land–sea carbon exports in the region.

Figure 3Modeled annual river loads of DIP (a), DSi (b), DIC (c), Alk (d), DOM (e) and POM (f).

Table 4Regional hotspots of riverine C loads (Tg C yr⁻¹) and DIP loads (10⁹ g P yr⁻¹) compared with regional estimates. Modeled DIP from our approach represents pre-industrial fluxes, whereas the DIP literature estimates are from present-day data and are strongly affected by anthropogenic perturbations.

¹ Araujo et al. (2014). ² Bird et al. (2008). ³ Tank et al. (2012). ⁴ Raymond et al. (2007). ⁵ Dittmar and Kattner (2003). ⁶ Le Fouest et al. (2013). ⁷ Li and Bush (2015). ⁸ Yoshimura et al. (2009). ⁹ Tao et al. (2010). ¹⁰ Seitzinger et al. (2010).

Download Print Version | Download XLSX

4.1 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 (Table 5). This implies that in REF, the amounts of P, N and Si added to maintain a plausible ocean biogeochemical state were similar to what is derived in our approach to estimate pre-industrial riverine exports. Despite slightly larger P and N inputs in RIV than in REF, RIV simulates lower global surface dissolved P and N concentrations, and lower global NPP. The coastal 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 water column. This is confirmed by the increased organic matter flux to the global coastal sediment (defined as areas with less than 250 m depths), with an increase from 0.18 Gt C yr⁻¹ (REF) to 0.25 Gt C yr⁻¹ (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⁻¹), 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 C inputs are on the other hand increased by almost 100 % in RIV in comparison to REF. The inorganic (366 Tg C yr⁻¹) and organic (237 Tg C yr⁻¹) C inputs in RIV show stronger agreement with the riverine inorganic (260–550 Tg C yr⁻¹) and organic C (270–350 Tg C yr⁻¹) 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; Regnier et al., 2013), signifying an improvement in the model's C inputs with respect to REF. These higher C inputs also result in a net long-term outgassing flux (231 Tg C yr⁻¹), 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₃ 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₂ and CO₂ outgassing in the long term (see Appendix C2).

Table 5Comparison of global oceanic inputs and global ocean biogeochemical state variables for REF and RIV. Additionally, we compare the modeled global mean surface DIP, DIN and DSi concentrations with World Ocean Atlas 2013 (WOA) surface layer means. OMZs: oxygen minimum zones.

Download Print Version | Download XLSX

The simulated global mean surface DIP and DIN concentrations are lower in RIV than in the observational data of the World Ocean Atlas 2013 (WOA) (0.439 and 0.480 µM DIP and 3.90 and 5.04 µM DIN respectively; WOA, 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 data set 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 data set 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 due to denitrification. Furthermore, the lower surface concentrations of DIP and DIN compared to those found in the WOA data set 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.

The DIP and DIN underestimation bias with respect to the WOA data set are also reflected in the spatial distributions of the surface concentrations (Fig. 4), in which the DIN concentrations are notably underestimated in most major basins. The spatial patterns of differences with regard to WOA data are similar for RIV and for REF (Appendix D, Figs. D1, D2, D3), suggesting that the physical ocean features are the dominant drivers of the nutrient distributions in the open ocean and of their bias (see also Supplement S5). Prominent spatial bias of both REF and RIV are underestimated surface DIP and DIN concentrations in the Southern Ocean, overestimated DSi concentrations in the Southern Ocean, and overestimated DIP concentrations in the tropical gyres.

Figure 4Surface DIP (a, d, g), DIN (b, e, h) and DSi (c, f, i) concentrations in WOA observations (a, b, c), RIV (d, e, f) and their differences RIV-WOA (g, h, i).

4.2 Riverine-induced NPP hotspots

The open ocean regions with the highest NPP rates in REF remain the most dominant areas of NPP in RIV (Fig. 5a, b), which also confirms the major importance of the ocean physics in dictating spatial patterns of the global NPP. In comparison 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 (Fig. 5c, d). This can be observed on the shallow shelves of the western tropical Atlantic, where major nutrient supplies are provided (see Sect. 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 high chlorophyll concentrations (Behrenfeld and Falkowski, 1997).

In the tropical Atlantic, tropical Pacific, North Pacific and Southern Ocean, the NPP is decreased in RIV with respect to REF, signaling that stabilizing ocean biogeochemical inventories with open ocean inputs, as was done in REF, also led to a slight artificial enhancement of the NPP in these regions. The decrease in the equatorial Pacific could partly be explained by the characteristics of South American river systems, which almost solely debouch into the Atlantic basin (Fig. 7). Although Southeast Asian rivers deliver substantial biogeochemical loads to the ocean, the export of the biogeochemical compounds from eastern Asian coastal systems to the open Pacific is likely inefficient, with modeled salinity profiles of the East China Sea suggesting little mixing with the open ocean (Supplement S4). Currents parallel to the coast could be a 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 the 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 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.

Figure 5Depth integrated annual NPP (a, b) and total annually accumulated organic matter concentration (c, d) in the surface layer in the RIV simulation and RIV – REF.

4.3 Riverine-induced CO₂ Outgassing

The addition of riverine C loads causes a simulated net oceanic CO₂ source of 231 Tg C yr⁻¹ to the atmosphere (Table 5). While the hotspots of the riverine-induced C outgassing are regions in proximity to major river mouths (Fig. 6b), a widespread, albeit weaker, outgassing signal can be observed in open ocean basins. The largest outgassing fluxes are found in the Atlantic and Indo-Pacific (31 % and 43 % of global outgassing flux respectively), which is due to large riverine C supplies from tropical Atlantic and Southeast Asian catchments (see Sect. 3.2). In the Southern Ocean, we observe an outgassing flux of 17 Tg yr⁻¹ when comparing RIV to REF, which is almost 10 % of the global riverine-caused outgassing. The riverine-induced outgassing is 113 Tg yr⁻¹ (49 % of the global outgassing flux) in the Southern Hemisphere, despite riverine C exports of the hemisphere contributing to only 199 Tg (33 %) of annual C loads to the ocean. This implies the presence of an oceanic interhemispheric transfer of riverine C from the Northern Hemisphere to the Southern Hemisphere (Fig. 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 a net outgassing of riverine C on the Laptev shelf and in the Beaufort Sea (2.2 and 2.3 Tg C yr⁻¹, respectively). The impacts of riverine C loads on the Arctic shelves can also be observed in the present-day coastal ocean pCO₂ data set of Laruelle et al. (2017), in which very high pCO₂ values are reported near river mouths (> 400 ppm).

Figure 6Annual pre-industrial air–sea CO₂ flux of (a) RIV, (b) RIV-REF, and (c) contributions of the Northern Hemisphere and Southern Hemisphere to riverine C loads and oceanic riverine C outgassing. A positive flux in (a) and (b) describes an outgassing flux from the ocean to the atmosphere, whereas a negative flux is a flux from the atmosphere to the ocean.

7.1 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 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 rivers could all be represented more realistically. We notably assumed a fixed riverine N:P mole ratio of 16:1 for all species and all catchments in this study. Riverine N:P ratios are suggested to exceed this ratio at the global scale (29:1 by Seitzinger et al., 2010 and 21:1 by Beusen et al., 2016), but denitrification in estuaries could compensate the excess of N with regard to a N:P ratio of 16:1 (see Supplement S1.3). Denitrification in estuaries is impossible to represent in the ocean model due to grid resolution constraints. Furthermore, we assume a similar spatial distribution of pre-industrial anthropogenic P inputs as for the present day. While the chosen spatial distribution agrees well with the distribution of other approaches (Supplement S1), we acknowledge a large degree of uncertainty both in the magnitudes of the pre-industrial anthropogenic inputs as well as in their distributions.

7.2 Dynamics of terrestrial organic matter in the ocean

The composition of terrestrial organic matter and its fate in the ocean have been strongly debated in the past. Recent 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 considered in this study, a differing C:P ratio to the one chosen in this study could affect the modeled 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 compounds.

7.3 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 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.