Implementing riverine biogeochemical inputs in ECCO-Darwin: a sensitivity analysis of terrestrial fluxes in a data-assimilative global ocean biogeochemistry model

Terrestrial sources of carbon and nutrients drive biogeochemical cycles in coastal regions and in the global ocean. Quantifying their impact on the spatiotemporal variability of the ocean carbon cycle is pivotal to understanding the distinctive characteristics of ocean basins dominated by riverine inflow. ECCO-Darwin is a data-constrained, global-ocean biogeochemistry model that has heretofore lacked lateral inputs of carbon and nutrients. The objective of this study is to add this new capability to ECCO-Darwin and to carry out a suite of sensitivity experiments in order to quantify the impact of these lateral fluxes on coastal- and open-ocean biogeochemistry. In this work, we use an optimized version of the data-assimilative global-ocean biogeochemistry ECCO-Darwin model to perform a sensitivity analysis of the ocean to lateral inputs of carbon and nutrients. We generate riverine inputs by combining daily point-source freshwater discharge from JRA55-do with the Global NEWS 2 watershed model, accounting for lateral inputs from 5171 watersheds worldwide. The addition of riverine inputs drives a small CO₂ outgassing (+0.02 Pg C yr⁻¹) due to compensating processes at regional scales. In basins dominated by carbon runoff, such as the Tropical Atlantic and Arctic Oceans, the addition of riverine inputs increases CO₂ outgassing (+13 % and +9 %, respectively). In contrast, runoff in nutrient-dominated Southeast Asia leads to increased CO₂ uptake (+9 %). This new riverine biogeochemical input capability will enable future ECCO-Darwin solutions to better capture key processes that occur along coastal margins in global oceans.

1 Introduction

Rivers transport carbon from land to the ocean as Dissolved Organic Carbon (DOC), Dissolved Inorganic Carbon (DIC), Particulate Organic Carbon (POC), and Particulate Inorganic Carbon (PIC), along with nutrients such as phosphorus, nitrogen, and silica, which are essential for phytoplankton growth. Terrestrial inorganic carbon and nutrients in streams originate from weathering of the lithosphere and the associated uptake of atmospheric CO₂, along with the remineralization of organic matter in streams and/or on land (Suchet and Probst, 1995; Battin et al., 2023).

Riverine carbon (0.7–1 Pg C yr⁻¹; Lacroix et al., 2021b; Resplandy et al., 2018; Gao et al., 2024; Liu et al., 2024) can be buried in coastal sediments, transported into the open ocean, and outgassed back to the atmosphere in the form of CO₂ (Liu et al., 2024; Regnier et al., 2022; Battin et al., 2023; Gao et al., 2024). This carbon is transferred to the atmosphere due to the saturation of surface-ocean waters by terrestrial DIC and the remineralization of terrestrial organic matter (Hartmann et al., 2009; Lacroix et al., 2020; Bertin et al., 2023) in shallow, well-mixed water columns. On continental shelves, the outgassing of CO₂ driven by the saturation of surface waters with terrestrial DIC or remineralized terrestrial organic carbon can also be compensated by the excess of alkalinity relative to DIC concentration (Cai, 2011; Louchard et al., 2021). In the absence of transformation in the coastal ocean, refractory riverine organic carbon can be transported offshore due to its slow turnover time (Hansell et al., 2004; Holmes et al., 2008; Kaiser et al., 2017). Concerning nutrients, their injection into the surface ocean can fertilize the growth of photosynthetic organisms in nutrient-limited regions. The subsequent primary production by photosynthetic organisms enhances CO₂ uptake by carbon fixation. Globally, lateral inputs increase ocean primary productivity and may contribute to an estimated coastal-ocean carbon sink from 0.2 to 0.7 Pg C yr⁻¹, which is roughly 10 % to 35 % of the global-ocean sink (Dai et al., 2022; Resplandy et al., 2024).

While monitoring global riverine inputs to the ocean is challenging due to the substantial financial/human effort, often in remote environments, land surface and watershed models can provide spatiotemporally-resolved lateral inputs at global scales (Mayorga et al., 2010; Krinner et al., 2005; Hagemann and Dümenil, 1997; Hagemann and Gates, 2003; Li et al., 2017; Bloom et al., 2020; Gao et al., 2023). Coupled with Global-ocean Biogeochemical Models (GOBMs), it is thus possible to quantify the response of the coastal- and open-ocean carbon cycle to lateral inputs (Aumont et al., 2001; Lacroix et al., 2021b; Mathis et al., 2022; Louchard et al., 2021; da Cunha and Buitenhuis, 2013; Le Fouest et al., 2013; Terhaar et al., 2019; Gao et al., 2023; Bertin et al., 2023; Manizza et al., 2019; Séférian et al., 2020). Here, we add the capability to represent lateral fluxes of carbon and nutrients in the ECCO-Darwin global-ocean biogeochemistry model and we examine the impact of these fluxes on the model’s sea-air CO₂ flux and Net Primary Production (NPP) state estimate to perform perturbation experiments attributed to lateral inputs of carbon and nutrients. ECCO-Darwin combines (i) property-conserving physics and circulation from the Estimating the Circulation and Climate of the Ocean (ECCO) project, (ii) the MIT Darwin Project’s marine ecology model, (iii) ocean carbon chemistry, and (iv) data assimilation tools developed by ECCO. The system provides global, data-constrained estimates of circulation, sea ice, ecology, and biogeochemistry, with demonstrated skill in reproducing variability in the carbon cycle (Carroll et al., 2020, 2022; Bertin et al., 2023).

In this study, we (1) add point-source lateral inputs of carbon and nutrients to ECCO-Darwin globally and (2) evaluate the model response of sea-air CO₂ flux and primary production to riverine inputs during 2000–2019. The sensitivity analysis described herein will allow for further understanding of the contribution of riverine inputs in future ECCO-Darwin solutions and ocean modeling studies that aim to represent processes occurring along coastal margins.

2 Methods

2.1 The ECCO-Darwin Ocean Biogeochemistry State Estimate

The ECCO-Darwin ocean biogeochemistry state estimate is extensively described in Brix et al. (2015), Manizza et al. (2019), and Carroll et al. (2020, 2022, 2024). For the ECCO-Darwin model presented in this study, ocean physics (circulation, temperature, salinity, and sea ice) are provided by a prerelease of the ECCO Version 4 release 5 (V4r5) global-ocean and sea-ice data synthesis. A detailed overview of ECCO V4 is available in Forget et al. (2015) while specific details pertaining to V4r5 are being made available in ECCO et al. (2024).

Horizontal discretization is based on a Lat-Lon-Cap-90 (LLC90) configuration of the MIT general circulation model (MITgcm; Marshall et al., 1997a, b). Nominal horizontal grid spacing is 1° but telescopes to ∼ 33 km meridionally near the Equator and to ∼ 55 km in the Arctic Ocean. The vertical discretization consists of 50 z-levels, ranging from 10 m thickness in the top 7 levels to 450 m at maximum depth of 6 km. ECCO V4 uses a third-order, direct-space-time tracer advection scheme in the horizontal and an implicit third-order upwind scheme in the vertical; a time step of 3600 s is used. Vertical mixing is parameterized using the Gaspar–Grégoris–Lefevre (GGL) mixing-layer turbulence closure and convective adjustment scheme (Gaspar et al., 1990). ECCO V4 assimilates physical observations via the adjoint method (Wunsch et al., 2009; Wunsch and Heimbach, 2013). Importantly, ECCO V4 is a property-conserving ocean reanalysis, that is, contrary to reanalyses that are based on sequential estimation methods, ECCO V4 satisfies model equations exactly for the complete period of optimization (1992–2020 for V4r5). This characteristic makes ECCO V4 uniquely well-suited for ocean ecology and biogeochemistry applications.

Daily river discharge in the present configuration is based on the Japanese 55-year atmospheric Reanalysis (JRA55) for driving ocean–sea-ice models (JRA55-do). JRA55-do river discharge is computed based on the Catchment-based Macro-scale Floodplain (CaMa-Flood) global river routing model and on adjusted runoff from the land component of JRA55 (Suzuki et al., 2018; Tsujino et al., 2018; Feng et al., 2021). JRA55-do point source freshwater runoff was added to ECCO V4r5 as a freshwater flux in the surface ocean (first vertical level) at the closest corresponding ECCO V4r5 grid cell along the coastal periphery. The freshwater flux was adjusted according to the difference in grid cell area between JRA55-do (0.25°×0.25°) and ECCO V4r5. A complete evaluation of ocean physics from ECCO V4r5 compared to observations can be found in the Supplement and in Feng et al. (2021).

ECCO V4r5 ocean physics were coupled online with the MIT Darwin Project ecosystem model described in Brix et al. (2015). The ecosystem model solves 39 prognostic variables, including carbon, nitrogen, phosphorus, iron, silica, oxygen, and alkalinity. The model simulates their respective cycle from inorganic pools to living/dead matter of plankton organisms and the subsequent remineralization, all driven by the ocean physics. The carbonate chemistry is solved with the method of Follows et al. (2006). Plankton species consist of five large-to-small functional phytoplankton types (diatoms, other large eukaryotes, Synechococcus, and low- and high-light adapted Prochlorococcus) and two zooplankton types. In the absence of lateral fluxes, carbon in ECCO-Darwin is removed from the ocean through a combination of biological, chemical, physical, and air–sea exchange processes. Phytoplankton uptake of DIC during photosynthesis reduces upper-ocean carbon and forms organic matter, some of which sinks out of the mixed layer as export production. Additional CO₂ drawdown occurs when surface waters are undersaturated relative to the atmosphere, leading to net air–sea CO₂ uptake. Carbonate chemistry processes, such as precipitation and dissolution, modify alkalinity and buffer the partitioning of carbon species, thereby influencing surface-ocean DIC concentrations. Finally, physical transport through upwelling, mixing, subduction, and advection transports both DIC and organic carbon through the water column. Nutrients are supplied by upwelling and vertical mixing, consumed by phytoplankton growth, regenerated during remineralization, and exported with sinking organic matter – collectively regulating the efficiency of carbon uptake and storage. In the water column, particulate matter (detritus, inorganic carbon, and living phytoplankton and zooplankton) sinks at prescribed velocities and is removed at the ocean bottom to limit the accumulation of particulates on the seafloor.

Atmospheric CO₂ partial pressure at sea level (ApCO₂) from the National Oceanic and Atmospheric Administration Marine Boundary Layer Reference product (Andrews et al., 2014) was used to drive sea-air CO₂ fluxes computed by the model according to Wanninkhof (2014). Atmospheric iron dust is deposited at the ocean surface based on the monthly climatology of Mahowald et al. (2009). ECCO-Darwin assimilates biogeochemical observations using a Green's Functions optimization approach (Menemenlis et al., 2005); the optimization methodology and associated data constraints are extensively described in Carroll et al. (2020). The ECCO-Darwin solution was previously published using an LLC270 ($\mathrm{1}/\mathrm{3}$°) ECCO solution (Zhang et al., 2018) and monthly climatological freshwater runoff forcing from Fekete et al. (2002). Here, we introduce a new 1°-version of ECCO-Darwin with daily point-source freshwater runoff from January 1992 to December 2019 (hereinafter our “Baseline” simulation) and also conduct a suite of perturbation experiments (Table 1) where we add various riverine biogeochemical input components to assess the primary productivity and carbon cycle response. Except for these changes, our simulations use the same initial conditions, parameter settings, and forcings as in Carroll et al. (2020). To account for biogeochemical spin-up in the perturbation runs, the following analysis was performed for the last 20 years of simulation, from January 2000 to December 2019 (Figs. S2–S9 in the Supplement).

Table 1 Sensitivity experiments and associated solutes: terrestrial dissolved organic carbon (t_DOC), dissolved inorganic carbon (t_DIC), dissolved inorganic nitrogen (t_DIN), dissolved organic nitrogen (t_DON), dissolved silica (t_DSi), dissolved inorganic phosphorus (t_DIP), dissolved organic phosphorus (t_DOP), dissolved inorganic iron (t_DFe), and dissolved organic iron (t_DOFe).

2.2 Baseline Evaluation

We compared simulated surface-ocean partial pressure in CO₂ (pCO₂) and sea-air CO₂ fluxes in Baseline with state-of-the-art products based on the Surface Ocean CO₂ Atlas (SOCAT; Bakker et al., 2016; Sabine et al., 2013). We used the monthly pCO₂ and sea-air CO₂ fluxes MPI-SOM-FFN v2023 (Landschützer et al., 2016; Jersild et al., 2023) and Copernicus Marine Environment Monitoring Service (CMEMS; Chau et al., 2022) climatologies computed from neural network-based clustering algorithms. In addition, we used the monthly atmospheric CO₂ inversion Jena Carboscope v2023 (Rödenbeck et al., 2013) based on high-precision measurements from the Gridded Fossil Emissions Dataset (GridFED; Jones et al., 2021) and SOCAT (Bakker et al., 2016; Sabine et al., 2013). These products were interpolated on the LLC90 grid from January 2000 to December 2019. Grid cells covered by sea-ice (concentration >0 %) were discarded from the model-data evaluation, based on the percentage of sea-ice cover simulated by ECCO-Darwin.

2.3 Biogeochemical River Discharge Product

In addition to the Baseline simulation, we conducted three sensitivity experiments (Table 1) where we added terrestrial DOC (t_DOC), DIC (t_DIC), total alkalinity (t_ALK), dissolved inorganic nitrogen (t_DIN), dissolved organic nitrogen (t_DON), and dissolved silica (t_DSi), dissolved inorganic phosphorus (t_DIP), dissolved organic phosphorus (t_DOP), dissolved inorganic iron (t_DFe), and dissolved organic iron (t_DOFe) henceforth referred to as riverine inputs in this study. Except for t_DIC, t_ALK, t_DFe and t_DOFe, riverine inputs are provided by the Global Nutrient Export from WaterSheds 2 (NEWS 2; Mayorga et al., 2010) model. The method for computing our daily point-source inputs, which is then used as forcing in ECCO-Darwin along the coastal periphery of the global ocean, is detailed below. Global NEWS 2 uses statistical and mechanistic relations at the watershed scale to compute annual-mean freshwater discharge and riverine inputs based on natural and anthropogenic sources, with 6292 individual watersheds delineated according to the global river systems dataset from Vörösmarty et al. (2000). Global NEWS 2 t_DIN was partitioned into nitrite (NO${{}_{\mathrm{2}}}^{-}$), nitrate (NO${}_{\mathrm{3}}^{-}$), and ammonium (NH${}_{\mathrm{4}}^{+}$), according to the mean fraction of each species concentration relative to the total DIN concentration from the GLObal RIver CHemistry Database (GLORICH; Hartmann et al., 2014). The NO${}_{\mathrm{2}}^{-}$ : DIN, NO${}_{\mathrm{3}}^{-}$ : DIN, and NH${}_{\mathrm{4}}^{+}$ : DIN ratios were estimated to be 0.02, 0.65, and 0.33, respectively. Inorganic phosphorus was partitioned into dissolved inorganic phosphorus (DIP) and iron-bound (Fe-P) pools using a fixed 1:3 DIP : Fe-P ratio based on pre-industrial export estimates (Compton et al., 2000). Fe was coupled to P at a $\mathrm{1}:\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ molar ratio, but the iron associated with the Fe-P oxide fraction was treated as non-bioavailable (Lacroix et al., 2020).

t_DIC inputs were computed using an empirical relation between freshwater discharge and gross CO₂ consumption from rock weathering, as described in Li et al. (2017, equation 9). CO₂ consumption by rock weathering over each Global NEWS 2 watershed was estimated based on the freshwater discharge and the basin-dominant lithology (Amiotte Suchet et al., 2003). t_ALK inputs were computed using an ALK : DIC ratio (0.98) based on the mean total ALK compared to DIC from GLORICH. The remineralization rate for terrestrial and marine DOC equals 1 over 100 d. We used Global NEWS 2 outputs for the year 2000 as representative of present-day carbon and nutrient inputs (Mayorga et al., 2010). Riverine inputs were compared against observations from literature and the Arctic Great Rivers Observatory (ArcticGRO) water-quality monitoring network in the Arctic region (Holmes et al., 2012; Tank et al., 2023) (Table S1 in the Supplement).

Global NEWS 2 river mouth locations were associated with JRA55-do grid points exhibiting the closest annual-mean freshwater discharge in 2000 within an euclidean distance of 5°. The top 100 largest rivers (by watershed extent) from Global NEWS 2 were imposed on JRA55-do grid points as a function of distance only. In total, 5171 river mouths were associated with JRA55-do grid points. For each discharge point, riverine input concentrations (g m⁻³) from the associated river were estimated by dividing the load by the annual volume of freshwater from Global NEWS 2; the concentration was then multiplied by the corresponding daily-mean freshwater flux from JRA55-do (m s⁻¹) to obtain a daily flux (g m⁻² s⁻¹). Riverine inputs were adjusted according to the grid-cell-area difference between JRA55-do and ECCO V4r5. Then, these biogeochemical inputs were added as point-source discharge along with riverine freshwater flux (Table 2 and Fig. S1). Due to overestimated t_DIC inputs in our Global NEWS 2-derived computation for the Amazon River, t_DIC inputs for this system were set to a more realistic, literature-mean of 2.54 Tmol C yr⁻¹ (da Cunha and Buitenhuis, 2013; Probst et al., 1994; Li et al., 2017) (for more details, see Appendix A). The outstandingly large Amazon watershed area (used for estimating rock weathering) and freshwater discharge compared to other basins drive a very high load when using equation 9 from Li et al. (2017). Riverine t_DOC, t_DIN, t_DON, t_DIP, t_DOP, and t_DSi inputs agree well with non-Global NEWS-based estimates in Table 2. t_DIC lateral inputs from rivers were estimated according to Amiotte Suchet et al. (2003), Mayorga et al. (2010), Li et al. (2017), resulting in t_DIC inputs of 381.81 Tg C yr⁻¹ to the ocean, which is in general agreement with recent studies in Table 2.

(Tian et al., 2023)(Li et al., 2017)(Liu et al., 2024)(Chen et al., 2025)(Tian et al., 2023)(Li et al., 2017)(Liu et al., 2024)(Sharples et al., 2017)(Ma et al., 2025)(Ma et al., 2025)(Turner et al., 2003)(Sharples et al., 2017)(Frings et al., 2016)(Turner et al., 2003)(Manizza et al., 2011)(Holmes et al., 2012)(Tank et al., 2012)(Sharples et al., 2017)(Holmes et al., 2012)(Holmes et al., 2012)(Sharples et al., 2017)(Sharples et al., 2017; Holmes et al., 2012)(Holmes et al., 2012)(Araujo et al., 2014)(da Cunha and Buitenhuis, 2013)(Araujo et al., 2014)(Sharples et al., 2017)(da Cunha and Buitenhuis, 2013)(Sharples et al., 2017)(da Cunha and Buitenhuis, 2013)(da Cunha and Buitenhuis, 2013)

Table 2 Riverine inputs and literature estimates from non-Global NEWS methods. NA – not available

2.4 Sensitivity Experiments and Analysis

Sensitivity experiments consisted of adding riverine inputs separately or together, along with freshwater runoff (Table 1). t_ALK was always added along with t_DIC in relevant experiments. Given that the previously optimized ECCO-Darwin solution did not include biogeochemical river discharge, the sensitivity experiments may contain some double-counting that will lead to deterioration of the model results relative to observed pCO₂ and sea-air CO₂ flux data products. Therefore, the analysis herein is restricted to examining the perturbation response rather than quantifying possible improvement or degradation of the simulation vs. observations. We analyzed monthly-mean model fields both in the coastal ocean (limits set by the furthest point from the coastline, either the 1000 m isobath or a distance of 300 km; 58×10⁶ km²) and open ocean (300×10⁶ km²) from 2000–2019. We also evaluated the sensitivity of ocean carbon cycling in three specific regions that receive large volumes of freshwater and biogeochemical inputs from major river systems (Lacroix et al., 2020): the Arctic Ocean (ARCT, 22×10⁶ km²), Tropical Atlantic (TROP-ATL, 77×10⁶ km²), and Southeast Asia (SE-ASIA, 62×10⁶ km²). Coastal and open-ocean boundaries are delineated by the black line shown in Fig. 1a.

3 Results

3.1 Baseline Evaluation

Overall, Baseline surface-ocean pCO₂ compares reasonably well with the Jena Carboscope, MPI-SOM-FFN, and Copernicus CMEMS data-based products (Fig. 1). The largest differences are concentrated along the coastal periphery and near large river mouths (i.e., Amazon, Paraná, Congo, Ganges, Yangtze, Amur), where Baseline underestimates surface-ocean pCO₂ (Fig. 1i). Additionally, the data-based products exhibited lower surface-ocean pCO₂ compared to Baseline (Fig. 1i) in the Arctic Ocean and near the periphery of Antarctica; regions where observations are highly limited in space and time.

Figure 1 Climatological global-ocean surface-ocean pCO₂ for (a) ECCO-Darwin Baseline, (b) Jena Carboscope, (d) MPI-SOM-FFN, (f) Copernicus CMEMS, and (h) mean of all data products. Panels (c), (e), (g), and (i) correspond to the difference between ECCO-Darwin Baseline and each data product. All fields shown are time means from January 2000 to December 2019. In (a), colored boundaries correspond to domains used for regional analysis of the Arctic Ocean (ARCT, violet line), Tropical Atlantic (TROP-ATL, red line), and Southeast Asia (SE-ASIA, green line). The black line delineates the coastal ocean from the open ocean, which is set by the furthest point from the coastline of either a 300 km distance or the 1000 m isobath. MPI-SOM-FFN, Jena Carboscope, and CMEMS Copernicus products were interpolated on the LLC90 grid.

Figure 2 Climatological global-ocean sea-air CO₂ flux for (a) ECCO-Darwin Baseline, (b) Jena Carboscope, (d) MPI-SOM-FFN, (f) Copernicus CMEMS, and (h) the mean of all data-based products. Panels (c), (e), (g), and (i) correspond to the difference between ECCO-Darwin Baseline and each product. Positive values represent CO₂ outgassing (red colors); negative values represent uptake (blue colors). All fields shown are time means from January 2000 to December 2019. In (a), colored boundary lines correspond to domains used for regional analysis of the Arctic Ocean (ARCT, violet line), Tropical Atlantic (TROP-ATL, red line), and Southeast Asia (SE-ASIA, green line). The black line delineates the coastal ocean from the open ocean, which is set by the furthest point from the coastline of either a 300 km distance or the 1000 m isobath. MPI-SOM-FFN, Jena Carboscope, and CMEMS Copernicus products were interpolated on the LLC90 grid.

Figure 2 shows a comparison of time-mean Baseline sea-air CO₂ flux (−2.58 Pg C yr⁻¹), Jena Carboscope v2023 (−2.11 Pg C yr⁻¹), MPI-SOM-FFN v2023 (−2.04 Pg C yr⁻¹), and Copernicus CMEMS (−1.97 Pg C yr⁻¹) products during 2000–2019. Compared to the product mean, Baseline sea-air CO₂ flux yields a stronger ocean CO₂ uptake (+0.5 Pg C yr⁻¹, Fig. 2i). Overall, the spatial distribution of source/sink patterns in the global ocean was well captured by Baseline. However, the products displayed stronger CO₂ outgassing in the North Pacific and Atlantic Oceans, the Arabian Sea, and the Southern Ocean (Fig. 2i).

3.2 Climatological Global Analysis

The addition of dissolved carbon and nutrients in ALL_run led to a small increase in CO₂ outgassing of 0.02 Pg C yr⁻¹ compared to Baseline, globally (Table 3 and Fig. 3a). The majority of CO₂ outgassing driven by riverine inputs (0.02 Pg C yr⁻¹) occurs in the coastal ocean (Table 3 and Figs. 3a and 4a). In ALL_run, the small net change in sea-air CO₂ flux results from compensation between the effects of riverine carbon and nutrients, as DC_run and NUT_run experiments result in elevated CO₂ outgassing and uptake, respectively (Table 3 and Fig. 4a). In DC_run, the increase in ocean carbon, and hence pCO₂ due to riverine inputs, reduces the ocean's capacity to take up atmospheric CO₂, resulting in a net CO₂ outgassing of 0.22 Pg C yr⁻¹ (Table 3 and Fig. 4a). In NUT_run, the increase of nutrients in the euphotic zone elevates phytoplankton productivity. The additional uptake of carbon by phytoplankton decreased surface-ocean DIC, resulting in an additional global-ocean CO₂ uptake of 0.20 Pg C yr⁻¹ (Table 3 and Fig. 4a).

Table 3 Sea-air CO₂ flux and Net Primary Production (NPP) for each experiment in the coastal, open, and global ocean. Positive values represent CO₂ outgassing; negative values are uptake.

While outgassing driven by carbon inputs was compensated by uptake due to nutrients in the open ocean, CO₂ uptake due to nutrients was 10 % lower than carbon-input-driven coastal outgassing, resulting in a global-ocean CO₂ uptake that was reduced by 0.02 Pg C yr⁻¹ (i.e., increased outgassing) in ALL_run compared to Baseline (Fig. 4a). Dissolved nutrient inputs in ALL_run resulted in a Net Primary Production (NPP) increase of 1 Pg C yr⁻¹ (+4 %) compared to Baseline (Table 3 and Figure 3b). The total increase of NPP in ALL_run from riverine inputs was stronger in the open ocean (0.62 Pg C yr⁻¹) compared to the coastal ocean (0.4 Pg C yr⁻¹) (Table 3 and Fig. 4b). However, the increase of NPP per surface area was larger in the coastal ocean (+7 g C m⁻² yr⁻¹, +7 %) compared to the open ocean (+2 g C m⁻² yr⁻¹, +1 %).

3.3 Climatological Regional Analysis

In Baseline, the CO₂ uptake in ARCT was roughly 0.21 Pg C yr⁻¹. When carbon and nutrient inputs are added in ALL_run, ARCT CO₂ uptake reduces by 0.02 Pg C yr⁻¹, with the majority of the response (75 %) in the coastal ocean (Table 4, Figs. 3a and 4a). In Baseline, ARCT NPP was 0.22 Pg C yr⁻¹, with a similar magnitude in the coastal and open ocean. Adding nutrient inputs into ARCT increased coastal NPP by 4 % (Fig. 4b).

Carbon and nutrient inputs resulted in a TROP-ATL CO₂ outgassing of 0.01 Pg C yr⁻¹ compared to Baseline (0.10 Pg C yr⁻¹). This imbalance results from CO₂ outgassing driven by dissolved carbon, which was 20 % larger than the uptake due to increased phytoplankton productivity from dissolved nutrients (Fig. 4b). In Baseline, NPP in TROP-ATL was 3.18 Pg C yr⁻¹. The increase in NPP driven by riverine nutrients occurs predominantly in the open ocean (∼ 65 %) compared to the coastal (∼ 35 %) zone (Figs. 3b and 4b).

SE-ASIA has a CO₂ uptake of 0.30 Pg C yr⁻¹ in Baseline, while CO₂ uptake increases by 0.02 Pg C yr⁻¹ in ALL_run (Table 4 and Fig. 3a). In the open ocean, the nutrient input-driven increase in NPP and associated CO₂ uptake is two times higher than carbon input-driven outgassing – leading to an overall imbalance and resulting in net CO₂ uptake in SE-ASIA (Fig. 4a). NPP in SE-ASIA without riverine inputs is 3.3 Pg C yr⁻¹. In ALL_run, NPP increases by 0.33 Pg C yr⁻¹ due to elevated nutrients in both the open and coastal ocean (Fig. 4b).

Figure 3 Global-ocean (a) sea-air CO₂ flux and (b) NPP driven by riverine inputs in ALL_run. In panel (a), positive values represent more CO₂ outgassing (red colors), and negative values represent more uptake (blue colors) compared to Baseline. Fields represent time-mean values from January 2000 to December 2019. Colored lines on maps show domains used for regional analysis. The black line delineates the coastal ocean from the open ocean.

Figure 4 Domain-integrated differences in (a) sea-air CO₂ flux and (b) NPP driven by exports in each sensitivity experiment. Differences were computed from time-mean fields from January 2000 to December 2019.

4 Discussion

4.1 ECCO-Darwin Baseline

Compared to state-of-the-art observation-based products, Baseline exhibits similar results to the version described in Carroll et al. (2020, 2022). Baseline depicts a time-mean global-ocean CO₂ uptake of 2.58 Pg C yr⁻¹ during 2000–2019. This is in relatively good agreement with MPI-SOM-FFN v2023 (−2.04 Pg C yr⁻¹; Landschützer et al., 2016; Jersild et al., 2023), Jena Carboscope v2023 (−2.11 Pg C yr⁻¹; Rödenbeck et al., 2013), and Copernicus CMEMS (−1.97 Pg C yr⁻¹; Chau et al., 2022) products over the same period. Lower Baseline surface-ocean pCO₂ and sea-air CO₂ fluxes compared to data-based products in the coastal periphery, especially near large river mouths, are driven by freshwater inputs only. In the absence of associated biogeochemistry, freshwater discharge dilutes chemical species in the coastal ocean, decreasing the salinity, the concentration of DIC, and the alkalinity in surface waters. This highlights the need to include coupled freshwater and biogeochemical discharge in GOBMs, as associated carbon and nutrients can compensate for the freshwater-only dilution effect. In the marginal ice zone of high latitudes, the data-based products depict lower surface-ocean pCO₂ and sea-air CO₂ fluxes compared to Baseline. As the data-products are primarily computed from statistical/mechanistic models based on the SOCAT database, the sparse observational coverage can be a source of error and uncertainty in these regions. We note that in regions such as the Antarctic Continental Shelf and the Arctic Ocean, which have extensive seasonal sea-ice cover, the SOCAT database coverage is limited (Bakker et al., 2016; Sabine et al., 2013).

Baseline captures similar spatial patterns of NPP compared to the model ensemble of the REgional Carbon Cycle Assessment and Processes Phase-2 (RECCAP-2) project that aims at constraining present-day ocean carbon from observation-based estimates, inverse models, and GOBMs (Doney et al., 2024) (Fig. S11). Many uncertainties remain regarding global-ocean NPP estimates from remote sensing (due to uncertainty in algorithms) and models (due to different conceptual model architectures). Overall, NPP in Baseline (24.5 Pg C yr⁻¹) lies in the lower bound of the wide range depicted by the RECCAP-2 model ensemble (25–57 Pg C yr⁻¹; Doney et al., 2024) and remote-sensing algorithms (43–68 Pg C yr⁻¹; Behrenfeld and Falkowski, 1997; Silsbe et al., 2016; Carr et al., 2006; Marra et al., 2003; Behrenfeld et al., 2005). This relatively low NPP results primarily from strong iron limitation in the High-Nutrient, Low-Chlorophyll (HNLC) regions in ECCO-Darwin (Carroll et al., 2020). The strong surface-ocean stratification and the weaker winter convection limit the replenishment of nutrients in the euphotic zone. Nevertheless, global-ocean NPP estimates will improve from enhanced space-time coverage of NPP measurements and associated key variables such as chlorophyll, light, nutrients, optical properties, and cell physiology (Bendtsen et al., 2023). An integration of environmental variables along with NPP measurements will greatly reduce models' spread and mismatch with synoptic in-situ observations. The implementation of a radiative transfer package (Dutkiewicz et al., 2019) in the next version of ECCO-Darwin, for which development is already underway, will permit the assimilation of direct ocean-color observations (remotely-sensed reflectance) and improve the model's estimate of global-ocean NPP.

4.2 Impact of Dissolved Carbon and Nutrient Inputs in ECCO-Darwin

We acknowledge that adding lateral inputs of freshwater, carbon, and nutrients in ECCO-Darwin Baseline can result in an additional source of spin-up and drift in the model simulations. As Baseline and sensitivity experiments are based on the same physical solution, the drift associated with the addition of freshwater is removed from our analysis; however, biogeochemical inputs may be an additional source of drift in the simulations presented in this study. The 28-year model period (1992–2019) does not allow the system to fully equilibrate with the addition of riverine inputs. However, time series of change in air–sea CO₂ flux and NPP with the addition of river carbon and nutrients (Figs. S2–S9) indicate that most regions approach quasi-equilibrium by the year 2000, consistent with the global response. In contrast, the change in air–sea CO₂ flux and NPP with the addition of river carbon and nutrients in the Arctic do not stabilize over the model period (Figs. S3 and S7). Regional variability in air–sea CO₂ flux responses can be interpreted through differences in coastal residence times, as in the Arctic, long residence times promote remineralization and outgassing of terrestrial organic matter while limiting nutrient-driven uptake due to light limitation (Liu et al., 2019; Lacroix et al., 2021a). These extended residence times also explain why the Arctic response does not stabilize within the 28-year experiment timescale (Figs. S3 and S7), in contrast to other regions where shorter residence times facilitate more rapid equilibration. Conversely, regions such as the Amazon plume display substantial CO₂ outgassing despite shorter residence times, but this is accompanied by elevated offshore transport, suggesting that riverine carbon inputs or remineralization rates may be overestimated in coastal systems where residence time is short. While the use of a Green's Functions-based optimization has been shown to reduce spin-up and drift in previous ECCO-Darwin solutions (Brix et al., 2015; Carroll et al., 2020), it will be necessary to optimize a new ECCO-Darwin solution that includes biogeochemical runoff to select the initial conditions and model parameters that will minimize model-data misfit (i.e., cost) and reduce spin-up drift – a focus of ongoing work. We note that the next version of ECCO-Darwin aims to include optimization controls of inputs ratio (DIC : ALK, NO${}_{\mathrm{2}}^{-}$ : DIN, NO${}_{\mathrm{3}}^{-}$ : DIN, and NH${}_{\mathrm{4}}^{+}$ : DIN), allowing us to optimize riverine inputs based on remotely-sensed and in-situ ocean observations.

Table 4 Change in sea-air CO₂ flux and NPP driven by riverine inputs. Positive values represent an increase in CO₂ outgassing or primary production; negative values represent an increase in CO₂ uptake or a decrease in primary production. NA – not available

^a Tivig et al. (2021); Cotrim da Cunha et al. (2007), ^b Manizza et al. (2011); Terhaar et al. (2021), ^c Effect of t_DOC only, ^d da Cunha and Buitenhuis (2013); Louchard et al. (2021), ^e Lower bound is for smaller domain in western TROP-ATL, ^f Tivig et al. (2021).

In this study, carbon inputs drive a CO₂ outgassing of 0.22 Pg C yr⁻¹, while nutrient inputs drive a CO₂ uptake of 0.20 Pg C yr⁻¹ from enhanced primary productivity, which primarily occurs in the coastal ocean. Combined, carbon and nutrient inputs in ALL_run are limited to an outgassing of 0.02 Pg C yr⁻¹ CO₂; lower than literature estimates (Table 4). In the simulation with riverine carbon only (DC_run), our estimate of +0.22 Pg C yr⁻¹ of air–sea CO₂ flux is lower than previous preindustrial-based estimates of +0.59 Pg C yr⁻¹ (Aumont et al., 2001) and +0.65 Pg C yr⁻¹ (Regnier et al., 2022), but within the same order of magnitude. When including both riverine carbon and nutrients (ALL_run), our model simulates a smaller increase in air–sea CO₂ flux (+0.02 Pg C yr⁻¹), alongside a positive NPP response (+1 Pg C yr⁻¹). This differs from the pre-industrial estimates of Lacroix et al. (2020), who found a comparable increase in air–sea CO₂ flux (+0.23 Pg C yr⁻¹) but a reduction in NPP (−1.78 Pg C yr⁻¹) due to stabilizing ocean biogeochemical inventories. Adding nutrient inputs increases global-ocean marine NPP by 1 Pg C yr⁻¹ compared to Baseline. The addition of t_DIN and t_DON also increased ocean NPP by 0.6 Pg C yr⁻¹ in the model described by Tivig et al. (2021) (Table 4). In our study, the increase in NPP per surface area driven by riverine inputs was stronger in the coastal ocean compared to the open ocean, relative to their respective surface areas. This is consistent with the recent study of Mathis et al. (2024), which demonstrates the role of increased nutrient inputs in driving stronger biological carbon fixation and, thus, an enhanced CO₂ sink in the coastal ocean during the last century. We note that our multi-decadal estimates do not reach equilibrium in the Arctic Ocean following the addition of riverine inputs (Figs. S2–S9) and do not have a realistic representation of blue carbon, bottom-sediment processes, and fine-resolution coastal ecosystems that drive the coastal-ocean sink and transformation of elements. Therefore, our results are not directly comparable to long-term and pre-industrial estimates of the ocean response to riverine inputs (Regnier et al., 2022; Resplandy et al., 2024).

Riverine inputs might be overlooked due to the lack of a more realistic representation of organic matter remineralization, allowing for the advection of excess dissolved carbon and nutrients into the open ocean. This may be due to our fixed DOC remineralization rate (100 d), which does not account for terrestrial-originating components with a faster degradation rate (labile to semi-labile), and the absence of a Land-to-Ocean Aquatic Continuum (LOAC) parameterization to account for estuarine and near-shore processes. For example, the strong CO₂ outgassing following the addition of riverine inputs on the Siberian Shelf in ALL_run may be driven by an excess of carbon reaching the ocean. Across the Arctic LOAC, permafrost DOC may be degraded and outgassed back to the atmosphere further upstream (river, estuary, river plume) compared to ECCO-Darwin, while our riverine inputs are directly injected into the coastal ocean (Spencer et al., 2015; Bertin et al., 2025). In TROP-ATL, DOC from the Amazon river is expected to be more stable in the coastal ocean (up to hundreds of years; Louchard et al., 2021). In SE-ASIA, excess inputs of nutrients reaching the Bay of Bengal or Sea of Japan may drive excess model perturbation in this region, as the model lacks a LOAC parameterization and especially representation of estuaries where nutrients can be consumed upstream by biological activity (Cai, 2011).

Assuming that carbon and/or nutrient inputs from each watershed are routed completely and instantaneously to the ocean is a source of model error, as losses and gains occur through the LOAC, especially in estuaries. Sharples et al. (2017) estimated that 25 % of global DIN inputs were removed on continental shelves through biological uptake, denitrification, and anaerobic oxidation. The absence of denitrification within estuaries (3–10 Tg N yr⁻¹) (Seitzinger et al., 2010) could alter N:P stoichiometry and downstream air–sea CO₂ fluxes. However, our results do not include sea-air CO₂ fluxes associated with these land-to-ocean components. Current GOBMs and Earth System Models (ESMs) used in IPCC Assessment Reports compute the amount of carbon introduced to coastal grid cells (i.e, lateral inputs) from reference watersheds or land-surface models that do not resolve the transport and transformation of carbon through the LOAC and, especially, estuaries and associated blue carbon pools (i.e., salt marshes and mangroves; Mayorga et al., 2010; Ciais et al., 2014; Lacroix et al., 2020; Ward et al., 2020). While coastal wetlands, estuaries, and continental shelves are a pivotal filter of carbon and biogeochemical elements, their action on reactive species has yet to be included in most GOBMs (Cai, 2011).

In Baseline, ARCT uptakes 213.9 Tg C yr⁻¹ of atmospheric CO₂; this may be an overestimate of the ARCT CO₂ sink, as recent estimates from modeling, atmospheric inversions, and pCO₂-based products range from 91–116 Tg C yr⁻¹ (Yasunaka et al., 2023) (Fig. S10); although observations are highly limited in this region. In ARCT, riverine inputs dominated by carbon reduce this CO₂ uptake by 20 Tg C yr⁻¹. In Terhaar et al. (2019), CO₂ outgassing increases by 90 % when riverine t_DOC was doubled. However, Terhaar et al. (2019) used an instantaneous remineralization rate for DOC, resulting in rapid outgassing in the coastal region compared to our results. In addition, nutrient inputs also contribute to the Arctic Ocean's carbon sink as they fertilize coastal waters. NPP in the Arctic Ocean increased by 4 % (+9 Tg C yr⁻¹) in ALL_run compared to Baseline. In Terhaar et al. (2019), the doubling of riverine nutrients (+2.3 Tg N yr⁻¹) leads to an 11 %-increase of NPP. Recent estimates by Terhaar et al. (2021) suggest that riverine nutrients support up to 15 % (+58 Tg C yr⁻¹) of marine NPP in the ARCT, in agreement with estimates by (Le Fouest et al., 2013, 2015) (Table 4). Therefore, biological CO₂ uptake driven by riverine nitrogen and its capacity to compensate CO₂ outgassing in ARCT might be underestimated in our study. We stress that the phytoplankton functional types in our global model are not representative of the specific Arctic Ocean ecology, and the lack of regionally-adjusted affinity for specific nutrients might hinder the model ecosystem response to riverine nutrients (Ardyna and Arrigo, 2020).

In Baseline, TROP-ATL is a source of CO₂ to the atmosphere (0.10 Pg C yr⁻¹), which agrees with both data-based products (Landschützer et al., 2016; Jersild et al., 2023; Rödenbeck, 2005) (0.04–0.08 Pg C yr⁻¹) and GOBM results (da Cunha and Buitenhuis, 2013; Louchard et al., 2021) (0.03–0.04 Pg C yr⁻¹). We note that previous studies show an input-driven increase in CO₂ uptake of 0.005 and 0.02 Pg C yr⁻¹ when adding biogeochemical runoff in TROP-ATL and western TROP-ATL, respectively (da Cunha and Buitenhuis, 2013; Louchard et al., 2021) (Table 4). However, in our simulations, the addition of riverine inputs in ALL_run enhanced the source of CO₂ to the atmosphere (+0.02 Pg C yr⁻¹) (Table 4). Contrary to the estimates of Louchard et al. (2021), which include physical effects associated with freshwater, such as enhanced upper-ocean stratification and gas solubility, our baseline simulation already includes these processes. Therefore, our set of experiments cannot isolate and quantify the impact of freshwater discharge on ocean biogeochemistry. Louchard et al. (2021) also included a regionally-adjusted plankton ecosystem, e.g., by including a nitrogen-fixing phytoplankton functional type, which increased the model's capability to resolve the biological pump and hence CO₂ uptake.

SE-ASIA is a sink of atmospheric CO₂ in Baseline (0.3 Pg C yr⁻¹). Combining sea-air CO₂ budgets for the different regions composed of SE-ASIA values from literature (East-Pacific, Indonesian seas, and North Indian Ocean (without including Oman and Somalian upwelling regions), we estimate an ocean carbon uptake of ∼ 0.2 Pg C yr⁻¹ for the entire SE-ASIA domain (Kartadikaria et al., 2015; De Verneil et al., 2023; Zhong et al., 2022; Hood et al., 2023). The net sea-air CO₂ exchange balance driven by riverine inputs in SE-ASIA results in a carbon uptake of 0.02 Pg C yr⁻¹ in ALL_run. Compared to ARCT and TROP-ATL, carbon uptake in SE-ASIA is enhanced by a large increase in marine NPP (+0.33 Pg C yr⁻¹, +9 %) driven by nutrient inputs. In Tivig et al. (2021), the simulated increase of NPP in response to riverine nitrogen was roughly 0.1 Pg C yr⁻¹ in Asia, with the strongest increase in the Yellow Sea, similar to our results (Table 4). Locally, adding riverine biogeochemical runoff also drives a source of CO₂ to the atmosphere, which is primarily limited to near river mouth locations in SE-ASIA. In the Yellow Sea and the Northern Bay of Bengal, close to the Yangtze and Ganges Rivers, the addition of riverine inputs at preindustrial levels in an ocean model also drove a CO₂ outgassing in Lacroix et al. (2020). Noticeably, in our simulations, the addition of carbon inputs switches the northern Bay of Bengal from a carbon uptake to a source, as suggested by Hood et al. (2023). Similarly, the addition of t_DOC in an ocean model of the Sunda Shelf Sea drives a CO₂ outgassing by 3.1 Tg C yr⁻¹ from 2013–2022 (Mayer et al., 2025). Most importantly, as nutrient inputs play a critical role in the SE-ASIA ocean carbon response, they need to be better constrained by a more extensive suite of observational data.

4.3 Recommendations for a More Realistic Representation of River-driven Carbon Cycling in ECCO-Darwin

This study presents a set of sensitivity experiments that quantify the contribution of riverine inputs in the ocean sea-air CO₂ flux and NPP; this was made possible following necessary and consequential simplifications that we elaborate in the following section. We also describe ongoing and future developments of ECCO-Darwin that will address these limitations and move toward a fully-optimized ECCO-Darwin solution that accounts for key processes along coastal margins.

t_ALK inputs were based on a global-mean, constant ALK : DIC ratio (0.98). We note that the GLORICH database used to compute this ALK : DIC ratio has relatively good coverage over the American continent; however, Eurasia and Africa remain underrepresented (Hartmann et al., 2014). As such, the ALK : DIC ratio can vary substantially over regional and time scales. The lack of this spatially-granular information in our simulated riverine inputs may misrepresent t_ALK inputs and the ALK-driven buffering capacity of simulated river plumes (Dubois et al., 2010; Tank et al., 2012; Mol et al., 2018; Ghosh et al., 2021; Gomez et al., 2023; Terhaar et al., 2019). While in estuaries, the absence of ALK relative to DIC leads to higher pCO₂ in upper-ocean waters and enhanced CO₂ outgassing in the coastal zone, rivers also result in an excess of ALK relative to DIC on continental shelves, which can reduce ocean pCO₂ through buffering and, thus, facilitate CO₂ uptake (Cai et al., 2010; Louchard et al., 2021).

In the present study, riverine particulate matter (1) rapidly sinks to the seafloor near river mouths, and (2) once at the seafloor, sinking particulates in the model are removed (at a rate equivalent to the sinking rate) to limit the unrealistic accumulation of particulates at depth. Remineralization of sinking particulates associated with riverine inputs and enhanced marine biomass could be an additional source of dissolved nutrients and carbon to the upper ocean through vertical mixing or upwelling mechanisms; ultimately affecting the sea-air CO₂ exchange depicted by the model in the coastal zone. In our current set-up, particulates from riverine-boosted production may be removed at the sediment-water interface too quickly, considering that most of the impact from riverine inputs occurs along the coast in shallow waters. Development to add a diagenetic sediment model in ECCO-Darwin is currently underway (RADI) to provide a more holistic representation of the global-ocean carbon sink (Sulpis et al., 2022).

Assuming that watershed-wide carbon and/or nutrient inputs are fully routed to the ocean is a misrepresentation, as losses and gains occur through the LOAC (Cai, 2011). Second, t_DOC is degraded in coastal waters at different rates depending on its origin and subsequent labile fraction (Holmes et al., 2008; Wickland et al., 2012; Shen et al., 2012; Lønborg et al., 2020). In the present study, in addition to not accounting for refractory and labile fractions of t_DOC, marine and terrestrial DOC are remineralized at the same rate (100 d). Overall, this could lead to unrealistic t_DOC remineralization in some regions and thus excess of either ocean CO₂ outgassing due to an excess of DIC or advection of organic matter to the open ocean; a limitation that also exists in other GOBMs due to undifferentiated remineralization rates. While recent modeling studies include separate pools of refractory and labile t_DOC with different remineralization rates at regional scales (Louchard et al., 2021; Gibson et al., 2022; Bertin et al., 2023), the nature of t_DOC needs to be better accounted for in GOBMs (such as in Aumont et al., 2001). For instance, the Amazon River – the largest global source of riverine t_DOC to the ocean – contributes to almost 50 % (+0.014 Pg C yr⁻¹) of the global-ocean CO₂ outgassing in our study. However, t_DOC from the Amazon River shows strong stability in the coastal ocean and is transported from the continental margin to the open ocean (Medeiros et al., 2015; Louchard et al., 2021). Increasing the refractory pool of Amazon t_DOC could, therefore, decrease CO₂ outgassing in our simulations. Nonetheless, the time scale of t_DOC remineralization remains difficult to constrain as observation-based estimates contain large variability in reported values (Holmes et al., 2008; Wickland et al., 2012; Shen et al., 2012).

5 Conclusion and Perspectives

In this study, we added the capability to represent lateral fluxes of carbon and nutrients in the data-constrained ECCO-Darwin global-ocean biogeochemistry model and we carried out a suite of sensitivity experiments in order to quantify the impact of these lateral fluxes on coastal- and open-ocean biogeochemistry. Globally, the role of present-day riverine inputs in ECCO-Darwin results in substantial, compensating regional responses in ocean carbon uptake and outgassing. In carbon-dominated margins, such as the Arctic and Tropical Atlantic Oceans, rivers drive a large source of CO₂ from the ocean to the atmosphere. In nutrient-dominated margins such as Southeast Asia, however, rivers drive a large ocean carbon uptake. While our experiments reveal clear regional responses, we identify limitations related to missing estuarine and benthic processing and incomplete equilibration over multi-decadal timescales in the Arctic Ocean. Our methodology combines Global NEWS 2 and JRA55-do to implement biogeochemical river discharge on top of point-source freshwater discharge globally, and at a daily frequency. These fields can be used (and are already being used) for many regional-to-global ocean model applications. Documenting such methodology is essential, given the lack of accurate representation of land-to-ocean and coastal processes in global ocean and Earth System Models (ESMs). This work is part of an open-science/open-source initiative available for everyone on the ECCO-Darwin GitHub repository (https://github.com/MITgcm-contrib/ecco_darwin/tree/master, last access: 20 January 2026). The quantification of the perturbation pertaining to the addition of terrestrial runoff in an ocean model over 20 years in the modern period is an interim, but significant step towards the development of new optimized ECCO-Darwin solutions that will integrate riverine inputs together with improved estuarine, sediment and benthic parameterizations.

Appendix A: Amazon River Runoff Set-up

As we computed riverine nutrient inputs from the combination of Global NEWS 2 loads with JRA55-DO runoff, Global NEWS 2 river concentrations must be co-located with JRA55-DO grid points exhibiting the closest annual discharge to avoid under- or overestimation of nutrient loads. In the case of the Amazon River, where freshwater and nutrient loads are extreme, we manually assigned the river mouth location from Global NEWS 2 to the corresponding JRA55-DO grid point. Additionally, when using equation in Li et al. (2017, equation 9), the DIC load from the Amazon river was overestimated and was therefore set to a mean literature value of 2.54 Tmol yr⁻¹ (da Cunha and Buitenhuis, 2013; Probst et al., 1994; Li et al., 2017).