Exploring Spatially Explicit Changes in Carbon Budgets of Global River Basins during the 20th Century.

Rivers play an important role in the global carbon (C) cycle. However, it remains unknown how long-term river C fluxes change because of climate, land-use, and other environmental changes. Here, we investigated the spatiotemporal variations in global freshwater C cycling in the 20th century using the mechanistic IMAGE-Dynamic Global Nutrient Model extended with the Dynamic In-Stream Chemistry Carbon module (DISC-CARBON) that couples river basin hydrology, environmental conditions, and C delivery with C flows from headwaters to mouths. The results show heterogeneous spatial distribution of dissolved inorganic carbon (DIC) concentrations in global inland waters with the lowest concentrations in the tropics and highest concentrations in the Arctic and semiarid and arid regions. Dissolved organic carbon (DOC) concentrations are less than 10 mg C/L in most global inland waters and are generally high in high-latitude basins. Increasing global C inputs, burial, and CO2 emissions reported in the literature are confirmed by DISC-CARBON. Global river C export to oceans has been stable around 0.9 Pg yr-1. The long-term changes and spatial patterns of concentrations and fluxes of different C forms in the global river network unfold the combined influence of the lithology, climate, and hydrology of river basins, terrestrial and biological C sources, in-stream C transformations, and human interferences such as damming.

Introduction

Rivers are an important component of the global carbon (C) cycle and have been identified as a significant source of carbon dioxide (CO2).[1] The estimates of global inland-water CO2 emissions range from 0.7 to 3.3 Pg C yr–1;[1−9] this large range not only implies uncertainty in global C budgets but also illustrates our limited understanding of governing factors.

C in freshwater originates from terrestrial (allochthonous or external) sources and in situ aquatic (autochthonous or within-system) production.[10] Allochthonous C is delivered to surface water as dissolved organic carbon (DOC) and particulate organic carbon (POC) from plant litter, surface runoff or leaching, and as dissolved inorganic carbon (DIC) produced during weathering or soil respiration.[11] After delivery to streams, rivers, lakes, or reservoirs, organic C can be metabolized to DIC, buried in sediment, or transported toward oceans.[1] The DIC delivered to or generated within the system is transported downstream or emitted to the atmosphere as CO2 since aquatic systems are predominantly supersaturated in CO2 relative to the atmosphere.[12−14] These early studies on CO2 partial pressure and effluxes have been an impetus to develop global assessments on the importance of freshwater systems in global C cycling.[1,8,15−17]

Most studies on C processing in streams, rivers, lakes, reservoirs, and floodplains have focused on local processes and cycling, including their sensitivity to perturbations, such as dam construction, eutrophication, land-use change, and climate change.[2,3,6,9,18−26] However, these available snapshot estimates often cover one single river, one specific year, or one specific C flux (mainly CO2 fluxes or river organic carbon export) and fail to describe the complete C budget. Moreover, these studies do not resolve how C cycling in aquatic systems has changed because of changes in hydrology, climate, and land use. Recently, Ran et al. (2021) reported a decrease in CO2 emissions from Chinese inland waters because of a combination of environmental factors,[27] and it is unknown whether such a change is specific to China or of global importance.

Modeling approaches are useful to describe long-term changes in the C cycle that occurred in the past and to generate projections into the future. Many existing river biogeochemistry models use regressions by lumping data at the river basin scale to quantify C export to coastal waters or CO2 emission to the atmosphere for specific years.[8,28−37] These lumped regression approaches are often based on a priori assumptions on the controlling factors or preselected datasets.[8,28,29,31−38] Such models lack spatiotemporal inputs and dynamic hydrological constraints, thus fail to reflect temporal changes and only improve our understanding of the underlying processes at highly aggregated levels (e.g., at river basin level or at one point in time). A few semi-mechanistic models that use some distributed inputs and in-stream processes have attempted to simulate the long-term changes in the riverine fluxes of certain C forms. For example, in the regional terrestrial ecosystem model (TEM), the annual river DOC export to oceans in Arctic regions is simulated as leaching of DOC from soil in spring based on a series of environmental factors, which causes huge uncertainties in estimates and inapplicability to other rivers.[39] The TRIPLEX-hydrological routing algorithm model (TRIPLEX-HYDRA) couples natural organic C inputs from soil and river hydrology and uses DOC observations from the literature to estimate the global riverine DOC export during 1951–2015, but the crucial in-stream C biogeochemical processes and all anthropogenic impacts including land-use change, wastewater discharge, dam construction, and so forth are not considered.[40] The Dynamic Land Ecosystem Model (DLEM) includes multiple C forms and in-stream processes and has been applied to some rivers in North America, but the potentially important processes of aquatic production, C burial, sediment dynamics, CO2 exchange with the atmosphere in large waterbodies, and the impacts of floodplains and construction of dams and reservoirs are not included.[41−43] The regional model Integrated Catchments Model for Carbon (INCA-C) is designed to simulate the processes of dissolved C, especially DOC, in boreal and temperate river basins, but sediment dynamics, POC dynamics, and anthropogenic impacts other than land cover are not considered.[44,45] The regional model ORCHILEAK simulates terrestrial C inputs, in-stream respiration and DOC decomposition, and transport of DOC and CO2 along the terrestrial-aquatic continuum of the Amazon basin, but other C forms and dynamic processes are not included.[46,47] The regional Model of Organic matter Removal and Export for Dissolved Organic Carbon (MORE-DOC) simulates riverine DOC processes in the Yangtze River for the period 1980–2015, but it relies on numerous observational data for calibration and does not include other C forms and relevant processes.[48] The regional process-based model National Integrated Catchment-based Ecohydrology (NICE)-BGC was recently applied to simulate the in-stream processing for DIC, DOC, and POC in global 153 river basins for the period 1980–2015 with a 1° × 1° spatial resolution, but the land-use change before the year 1992, impacts of dams and reservoirs, and processes relevant to sediment and C burial were not considered.[49,50] Since these process-based models do not fully consider multiple C forms and their associated in-stream dynamic processes, including C production, consumption, transformation, lateral transport, and interaction at the water-sediment and water-air interfaces, they have limited capabilities of hindcasting the historical spatially-explicit riverine fluxes and concentrations of multiple C forms on the global scale.

To describe the changing riverine C cycling resulting from the long-term interactions between land-use changes, interventions in hydrology (dam construction, reservoirs, and water extraction), and wastewater discharge, we need a spatially explicit integrated model that describes biogeochemical processes coupled to hydrology and the terrestrial C cycle. In this study, we implement freshwater C cycling in the process-based Dynamic In-Stream Chemistry module (DISC) (i.e., DISC-CARBON), which is part of the Integrated Model to Assess the Global Environment (IMAGE[51])-Dynamic Global Nutrient Model (IMAGE-DGNM[52]). This new model describes the spatial and temporal variability of C concentrations, transformations and fluxes based on the river basin hydrology, environmental conditions, and C delivery from headwaters to mouths. The model calculates pCO2, CO2 emissions, organic carbon burial, and export of DIC, POC, and DOC resulting from the balance of inputs, transfers, and transformations, and observational data are only used for performance assessment and not for calibration. We employ DISC-CARBON to investigate long-term changes in concentrations and fluxes for DIC, DOC, and POC in the world’s river network over the 20th century, and, via a sensitivity analysis, identify the major drivers of C fluxes (export, burial, and emission) for five major river basins.

Model and Data Used

DGNM Framework

The DISC-CARBON module, part of the IMAGE-DGNM framework (Figure a), is an extension of the recently published DISC module[52] with a description of the riverine C cycle. IMAGE-DGNM builds on the IMAGE-GNM[53] that uses the spiraling approach for describing global annual in-stream nitrogen (N) and phosphorus (P) retention for long time series (20th century). With defined subannual variations and speciation of C and nutrient delivery to river networks,[54] and the integration of the new process-based DISC module that replaces the spiraling approach, the IMAGE-DGNM framework allows for the global simulation of transfers of multiple nutrient and C forms from land to coastal waters. IMAGE-DGNM provides a long-term global perspective on C and nutrient accumulation and consumption processes in landscapes within river basins. Apart from the subannual temporal scale and the representation of different nutrient and C forms, the description of remobilization processes of the accumulated matter is a major innovation at this scale of analysis (with more details in Text SI1 in the Supporting Information).

Figure 1

(a) Scheme of the IMAGE-DGNM framework including the DISC-CARBON module for the in-stream biogeochemical C transformation processes and (b) scheme of C sources, forms, and biogeochemical transformations in all simulated waterbodies in the DISC-CARBON module. The formulas for each transformation process are listed in Table SI2.

In the current version of the framework, input and output time steps range from monthly to yearly. In this study, we used an annual temporal resolution at the global scale and compared results for annual and monthly resolution for the Rhine River basin. The spatial resolution of 0.5-by-0.5 degree matches that of the PCRaster Global Water Balance (PCR-GLOBWB) hydrology model, which is part of the IMAGE-DGNM framework (Figure a). PCR-GLOBWB dynamically simulates the volumes, surface areas, and discharges of the different waterbodies of river networks, including lakes, reservoirs, and high-order (≥ 6) streams.[55,56] Floodplains are assumed to exchange water in Strahler orders ≥6 and have a flow velocity of 10% of that in the mainstream. The discharge and characteristics of smaller streams are estimated for each 0.5 by 0.5-degree continental grid cell based on the runoff and the properties of high-order streams using the parameterization proposed by Wollheim et al.[57]

Within the IMAGE-DGNM framework, the IMAGE model provides land cover data to PCR-GLOBWB and C delivery fluxes to DISC-CARBON. Climate data from ERA-40 re-analysis[58] are used in PCR-GLOBWB for computing the water balance, runoff, and discharge for each year. C delivery includes POC, DOC, DIC (the sum of dissolved CO2, HCO3–, and CO32–) and alkalinity (ALK) from wastewater, surface runoff, weathering, eroded soil material, and litterfall from vegetation in flooded areas (Figure a). The calculation of each of these input fluxes is described in Table SI1.

DISC-CARBON

After delivery of C to streams and rivers, the DISC-CARBON model calculates in-stream biogeochemistry, uptake by pelagic and benthic algae (ALG), burial, mineralization, CO2 emission, and transport for all waterbodies from upstream grid cells down to the river mouth for each river basin at the global scale (Figure b). Particulate inorganic carbon (PIC), primarily calcium carbonate, is ignored in DISC-CARBON for simplification, considering that PIC mainly originates from detrital old carbonates,[59] that PIC transport flux is relatively limited compared with those of other C species,[16,25,60] and that its relocation within river basins does generally not affect the short-term biogeochemical C fluxes.[61] DISC-CARBON calculates the concentration of C species i for each time step in each waterbody type (i.e., rivers of different stream orders, lakes, reservoirs, or floodplains) of every cell as an effect of biogeochemical (bgc) interactions between C species and as a result of hydrological (hyd) transport (between cells or waterbodies within the same grid cell) as follows:

if species i is dissolved or suspended in the water column:

if species i is incorporated or attached to the bed surface:

Hydrological advection of any dissolved or particulate C species i in the water column, being DIC, DOC, POC, and ALG, is calculated as follows:where L is the upstream load (Mmol yr–1) of species i, Q is the water discharge (km3 yr–1), and [C] is the concentration of C species i (Mmol km–3). POC, when sedimented, is not transported until it is resuspended. Benthic algae are assumed to be attached to the streambed and are not transported downstream.

Model equations for C dynamics are listed in Table SI2, and model constants and parameters are listed in Tables SI3 and SI4, respectively. Subsequently, we briefly discuss the simulation of light limitation, primary production, and DIC and DOC dynamics.

The biogeochemical transformations among DOC, POC, and DIC depend on hydrology, temperature, and radiation. Light limitation for primary production is calculated using a spatiotemporal distribution of solar radiation reaching the surface of the waterbody and water turbidity, caused by particulate matter, which controls light penetration through the water column (eqs 1–7 in Table SI2). Primary production (eqs 8–19 in Table SI2) depends on the biomass of the producers, their growth rates, temperature, and light and DIC availability. Similarly, respiration and excretion (biomass to DOC) are modeled as a fraction of primary producer biomass and depend on temperature.

POC dynamics is affected by delivery fluxes from erosion and litterfall, sedimentation and resuspension, primary production, and mineralization. The mineralization of terrestrial organic matter with structural carbohydrates and lignins is slower than that of aquatic organic matter, which is rich in N and P.[62] The role of associations with protective mineral surfaces is ignored (see previous studies[63−66]). However, to account for the diversity in POC reactivity,[67,68] DISC-CARBON distinguishes (1) allochthonous, terrestrial POC (with slow mineralization) and (2) aquatic, autochthonous POC (with fast mineralization) (Table SI2 eqs 27–30).

The dynamics of DIC is described by eqs 20–25 in Table SI2. External input of DIC in the DISC-CARBON module originates from weathering (Table SI1), and DIC is produced in-stream through the mineralization of organic C forms and respiration of living biomass. DIC consumption occurs through primary production. Finally, DIC is added to or removed from the waterbody through CO2 exchange with the atmosphere. ALK is generated by weathering of sediments and rocks and delivered to streams. ALK production and consumption by primary production, respiration, nitrification, calcium carbonate precipitation, and dissolution within the stream network[69] are assumed to be negligible compared to the ALK inputs from weathering; the spatial and temporal resolution of the model does not allow for calculating hyporheic ALK production and removal.[70] Although ALK is the sum of excess bases in solution in natural environments, carbonate alkalinity (i.e., the bicarbonate and carbonate ions) tends to make up most alkalinity. In DISC-CARBON, ALK delivered to surface water is combined with DIC to calculate pCO2 and pH.

Surface runoff and wastewater constitute external DOC sources (Table SI1). In-stream DOC production occurs through excretion by pelagic and benthic algae, and DOC consumption occurs through mineralization. DOC dynamics is described by eqs 26–46 in Table SI2.

Sensitivity Analysis

For a selection of the world’s largest rivers (the Amazon, Lena, Mississippi, Nile, and Yangtze), we calculated the sensitivity of the DISC-CARBON modeled average CO2 emissions, C burial, and C export [DOC, POC, DIC, ALG, and their sum (total carbon, TC)] to the variation of 52 parameters, 8 environmental, and 8 C delivery forcings. We used the Latin hypercube sampling method[71] and carried out 750 runs. Model sensitivity is expressed as the standardized regression coefficient (SRC). More details on this approach are provided in Text SI2.

Model Performance Assessment Strategy

The model performance is evaluated by comparing simulated concentrations of DIC, DOC, and TC with measurement-based annual and 0.5-by-0.5-degree grid averages for a range of global rivers calculated from GloRiSe,[72] GLORICH,[73] and GEMS-GLORI databases.[74] Available data have been selected for stations and years with at least six observations per year to represent the annual mean, while these measurements may include different sampling and analytical methods (which are not always recorded in the databases).

Results and Discussion

Model Performance

Model Estimates and Observations

The comparison of the simulated concentrations of DIC, DOC, and TC with the available observational data in river basins since the 1940s shows that the concentrations of different C forms simulated by DISC-CARBON are in fair agreement with observations (Figure a–c). Most simulated DIC (89%) and DOC (90%) concentrations are within a factor of 2 of the observations at various stations in a range of global river basins, irrespective of the basin size (Figure a,b,e and SI1a). The simulated TC concentrations agree with the measurements for major rivers; for small rivers, 83% of the simulated TC concentrations are within a factor of 2 of those measured, with a slight general overestimation (Figure c and SI1b). For data from river mouths, the simulated DIC (R2 = 0.519, p < 0.001) and DOC concentrations (R2 = 0.562, p < 0.001) show an even better agreement with observations, with over 97% of both within a factor of 2 (Figure SI2a, b). DISC-CARBON predictions of POC river export are consistent with observations for large global rivers (R2 = 0.462, p < 0.001) (Figure d). Despite the limited available measurements of TC (i.e., all C forms measured simultaneously at the same location) at global river mouths compared with those of DIC and DOC, 97% of the simulated TC concentrations are within a factor of 2 of the measurements (Figure SI2c). Freshwaters are heavily under-sampled and the lack of representativeness of the few TC measurements available may also contribute to the minor mismatch. DISC-CARBON results are closer to the observations for river mouths than for upstream stations because the simulation results at river mouths integrate the effect of processes over the entire river basins: the quality of the input data is better at the basin than the single-grid-cell scale. For upstream stations, the spatial input data of DISC-CARBON can be uncertain because of the coarse resolution, and available observational data may not be representative because of incomplete temporal coverage within each year.

Figure 2

Validation of DISC-CARBON against the observed concentrations of (a) DIC, (b) DOC, and (c) TC (DIC + DOC + POC) at different stations with at least six measurements within the year considered for numerous global rivers of various sizes; when more than one station occurs within a grid cell, the mean of their annual average concentrations is used for comparison; (d) validation of DISC-CARBON simulated POC export to the coastal oceans against observation data from the late 1980s to early 1990s; (e) fraction of observations plotted against the ratio of prediction: observation (relative error) for DIC; and (f) comparison of the difference between the predicted and observed DIC concentrations (in μmol/L) with the mean of the predicted and observed values according to Bland and Altman (1986).[75] Similar figures for DOC and TC to (e) and (f) for DIC are in Figures SI1a–d. DIC data covering 1942–2000, DOC data covering 1973–2000, and TC data (for Rhine and Weser Rivers) covering 1978–1998 are from GloRiSe[72] and GLORICH;[73] POC observation data are from GEMS-GLORI.[74] Rivers are sorted based on their catchment areas.

We also tested the model performance with the Bland–Altman approach,[75] that is, examining the difference between the observation and prediction (residual) with the mean of the predicted and observed values. The results show that the predictions agree with the observations (Figure e, SI1a, b), and there is no systematic error for DIC (Figure f), DOC (Figurer SI1c), and TC (Figure SI1d).

This fair agreement between model predictions and independent observations gives confidence to the overall approach because DISC-CARBON is a mechanistic model and not a regression model aiming to reduce differences between observations and estimates.

Sensitivity Analysis

Since our model is based on mass conservation, the modeled global C fluxes are internally consistent, and temporal changes are the combined result of the changes in hydrology, climate, and C inputs from the land. However, the uncertainties in the fluxes can be as large as their temporal changes, which calls for an analysis of the sensitivity of modeled fluxes to the variation of model parameters. The results of this sensitivity analysis for those parameters with a significant and important influence on the simulated CO2 efflux, C burial, and export of different C forms are presented in Table SI5 for the Amazon (a), Lena (b), Mississippi (c), Nile (d), and Yangtze (e) Rivers. These rivers represent a range in climate, hydrology, and human activity (land use, dam construction, and so forth). In this way, not only the variations of parameters driving changes in C cycling but also their differences among river basins are examined. We present values of the SRC, which quantify the sensitivity of the model to parameters (see Tables SI5). Results for 32 input parameters are shown that have a significant and important effect on one of the output variables listed in Table SI5 in any of the five rivers analyzed.

For the five rivers analyzed, the role of discharge (Q) is important for almost all C fluxes as it influences the flow velocity and time available for transformations, with a commonly negative influence of increasing discharge on CO2 efflux and C burial and a positive influence (although limited in some cases) on the export of the various C forms (because there is less burial and less CO2 efflux), which is consistent with the findings in other studies.[36,44,46,48,76,77] An increase in temperature positively influences CO2 efflux (SRC of 0.28 for the Amazon to 0.39 in the Yangtze and 0.65 in the Nile) and lowers C burial (not important in the Lena because of low temperatures and low production rates, and SRC values of −0.43 to −0.83 in the other rivers); it also negatively impacts the export of primarily DOC (high negative values of −0.43 for the Lena and −0.73 for the Amazon) and POC (−0.22 for the Nile to −0.73 for the Amazon). The parameters related to phytoplankton growth (e.g., solar radiation and maximum rates of phytoplankton growth or mortality) are important in most of the rivers for the export of POC and algal biomass (SRC values of 0.22–0.30, Tables SI5a–e), but they also exert an influence on DOC export, C burial, and CO2 emission in rivers heavily controlled by dams and reservoirs, such as the Mississippi[17,78] and Yangtze.[27,79,80] Maximum algal growth rates influence CO2 exchange negatively (SRC of −0.36 for the Mississippi and Yangtze) and burial positively (SRC of 0.48 for the Mississippi and 0.43 for the Yangtze); the sensitivity to the variation of parameter values for algal growth is much smaller in the other rivers considered (Tables SI5c and e).

Terrestrial organic inputs generally have a small positive effect on the POC export (highest SRC values of 0.41–0.43 for the Mississippi, Nile, and Yangtze, Tables SI5c,d,e). Differences in geomorphology play a role in the importance of autochthonous phytoplankton growth. For the Amazon River and its tributaries, floodplains are major components[81−83] and have large standing stocks of biomass and limited agricultural land use, which cause large C inputs to the water during flooding periods in forested areas (POC input from litterfall, with a high SRC of 0.85 for CO2 emission). In contrast, the Mississippi and Yangtze Rivers have extensive agricultural areas,[84,85] and floodplain processes are less important (SRC of 0.39 and 0.2 for CO2, respectively, Tables SI5c and e) than those in the Amazon. With its extensive forested area in upstream areas, the Nile basin also has a high SRC of 0.58 for CO2 emission because of variation in litterfall (Table SI5d). Finally, alkalinity inputs from groundwater have a significant and important influence on DIC export in all analyzed rivers (with SRC values ranging from 0.77 for the Amazon to 1.0 for the Rhine), and in most cases, there is also a large influence on the TC export, especially where DIC is the major component of TC in the rivers.

Although the small parameter ranges used to calculate the model sensitivity are by no means uncertainty ranges, the approach allows for calculating the ranges covering 95% of the outcomes. We illustrate the model behavior with the estimates and their ranges for the five rivers in Table . The results indicate that the ranges vary among rivers and also among output variables. The CO2 emissions vary less in the Amazon River than in the other rivers, the variation in C burial is the largest for the Yangtze River and that of C export is the largest in the Amazon and Yangtze Rivers. However, our analysis demonstrates that relatively narrow ranges of 5% for all inputs and parameters except temperature (Text SI2) can lead to large variations in the model output.

Table 1

Average Yearly CO2 Emission, TC Burial, and TC Export by the Amazon, Lena, Mississippi, Nile, and Yangtze Rivers over the Period 1995–2000 Simulated by DISC-CARBON with the Ranges Covering 95% of the Outcomes (between Brackets in %) Based on the Assumed Ranges in Input Parameters Listed in Text SI2

rivers	CO2 emission	TC burial	TC export
	Tg C yr–1
Amazon	975.2 (±5%)	33.2 (±15%)	89.9 (±13%)
Lena	2.6 (±11%)	0.03 (±13%)	4.9 (±4%)
Mississippi	20.7 (±11%)	18.1 (±8%)	19.1 (±5%)
Nile	1.6 (±11%)	0.1 (±11%)	2.7 (±5%)
Yangtze	2.6 (±16%)	0.1 (±21%)	13.4 (±13%)

DIC and DOC Concentrations in Global Inland Waters

After the above performance evaluation, we simulated C cycling in global river basins for the period 1900–2000 at a yearly time step. The results of dissolved C concentrations in global inland waters with a 0.5° × 0.5° spatial resolution are shown in Figure and SI3.

Figure 4

(a) Aggregated C burial, export, and emission (i.e., CO2 emission) in global river basins for the period 1900–2000 and (b) export of POC, DOC, DIC, and TC to world oceans in 1950 and 2000. C input is the sum of burial, emission, and export.

DIC Concentration Patterns

In 2000, simulated DIC concentrations in global inland waters show a wide range with distinct regions having high or low concentrations (Figure a). The lowest simulated DIC concentrations are in the equatorial region between 10°N and 10°S, with a range of 10–30 mg C/L in the island area of Oceania and southern Asia and even lower levels (minimum values <10 mg C/L) in the continental areas of Africa and South America. These results agree with the observed mean DIC concentrations of 8 mg C/L in the Amazon River and 3 mg C/L in the Congo River.[59] Huang et al. (2012) reported a mean DIC concentration of 8 mg C/L for 175 equatorial rivers;[25] our results are also consistent with the reported mean DIC concentrations of 5 mg C/L in Africa and America, 13 mg C/L in Asia, and 21 mg C/L in Oceania. The slightly higher DIC concentration in the equatorial region of Asia than in South America is mainly because of the more widespread presence of carbonate rock in Asia.[34,86]

Figure 3

Concentrations of DIC (a) and DOC (b) in global inland waters simulated by DISC-CARBON for the year 2000. Gray colors indicate grid cells with precipitation excess lower than 3 mm per year. Maps showing the concentration difference between 1950 and 2000 are presented in Figure SI3.

The highest predicted DIC concentrations (generally over 40 mg C/L) in global inland waters are mostly in the Arctic river basins (in eastern Asia, Europe, and North America), which are regions with carbonate and volcanic rocks and regions with calcareous soils (mainly loess).[24,25,86−88] The simulated spatial distribution of DIC concentration in global inland waters thus seems closely related to lithology, which is consistent with the early study by Meybeck et al.[89] and with the consensus that DIC in rivers primarily originates from rock and sediment weathering.[90,91] Therefore, driven by weathering of carbonate (in North America) and silicate rocks (in Siberian watersheds),[88] the high DIC loadings of Artic rivers may be an important source of DIC to the Arctic Ocean. The simulated high DIC concentrations in Arctic rivers also agree with observations.[92] Simulated DIC concentrations are also high in rivers in semiarid and arid regions because of the high evaporation and low water discharge,[93] including the Niger and Nile River basins (in northern Africa), Limpopo and Orange River basins (in southern Africa), Darling River basin (in Oceania), and Indus, Ural, Volga, Tigris, and Euphrates River basins (in eastern Europe and western Asia).

DOC Concentration Patterns

In 2000, simulated DOC concentrations in inland waters are usually lower than 10 mg C/L (Figure b), which is consistent with the mean DOC concentration in global rivers of 6 mg C/L[25] and reports of 6 mg C/L in Oubangui River, 5 mg C/L in Mpoko River, 10 mg C/L in Ngoko-Sangha River, and 11 mg C/L in Congo-Zaïre River during 1990–1996.[25] The DOC concentrations tend to be higher at high latitudes, especially in the Northern Hemisphere, and the highest DOC concentrations are in Arctic river basins. This agrees with the reported high DOC concentrations for the Yenisey (2–13 mg C/L), Ob’ (4–17 mg C/L), Lena (3–24 mg C/L), Yukon (3–16 mg C/L), Porcupine (2–12 mg C/L), and Kolyma (3–18 mg C/L) Rivers.[87,94,95] The latitudinal distribution of river DOC concentrations is mainly due to the spatial heterogeneity of soil C inputs.[94,96] High-latitude soils and peatlands account for about half of the global soil C stock, much of which is in the Arctic watersheds that extend as far south as 45°N in the Eurasian continent.[94,96] With substantial inputs from the large soil C stocks,[23,87,94−96] the DOC concentrations in Arctic rivers are among the highest. The DOC concentrations for tropical rivers are low because high DOC inputs (from vegetation and soils) to rivers are balanced by high DOC decomposition rates at high temperatures.

This consistency between model predictions and observations for DOC and DIC at the various sites discussed above indicates that DISC-CARBON sufficiently incorporates the combined influence of the lithology, climate and hydrology of basins, terrestrial and biological sources, and in-stream transformations on the distributions of DIC and DOC in global rivers.

Global River C Budget and Export

To evaluate the simulations of global C fluxes in inland waters and export to coastal waters, it is important not only to understand the governing factors at the global scale and how well they are represented in our models but also to assess how simulated global C fluxes compare with data-based estimates. No direct observations are available to validate the C input from terrestrial ecosystems, and there are a few observations for C burial and CO2 emissions. Moreover, these observations are only available for the last part of the twentieth century. We, therefore, evaluate the flux estimates, their changes, and uncertainties by comparison with estimates based on data, models, and other approaches from the literature.

The literature inventory shows that the temporal changes simulated by DISC-CARBON (Figure a) are generally consistent with the trends from earlier data-based and model-based assessments (Table ). A direct comparison of our model predictions with these assessments is not straightforward because in these studies, (1) “static” C fluxes are based on merging data covering long periods, (2) the data usually did not cover all global inland waters, and (3) the whole C budget was rarely fully quantified with a consistent approach, as done in our study. Literature estimates of C inputs from terrestrial ecosystems are basically the resultant of the overall C budget and are consequently subject to substantial uncertainty. Our estimate of the C input to global inland waters for the year 2000 (Figure a) is in the middle of the range of previous estimates for the 2000s of 1.9 Pg C yr–1 by Cole et al.[1] to 5.1 Pg C yr–1 by Sawakuchi et al.[6] and Drake et al.,[11] and is very close to 2.7 Pg C yr–1 estimated by Battin et al.[5] and Regnier et al.[9] and 2.9 Pg C yr–1 estimated by Tranvik et al.[2] (Table ). These estimates have been derived in multiple ways and are generally based on the extrapolation of independent C burial, CO2 emission, and river C transport data, unlike our integrated model results. Our estimated 3.0 Pg C yr–1 of C input to inland waters for the year 1900 is higher than 1.7 Pg C yr–1 for the preindustrial era estimated by Regnier et al.[9] This difference is likely due to multiple factors: for example, we included, while they ignored, the decline of natural inputs because of land-use changes over the past century (Figure SI4).

Table 2

Comparison with Existing Estimates of Global Freshwater Carbon Fluxes

study	period	TC export	CO2 emission	TC burial	TC delivery from land
		Pg C yr–1
Mulholland and Elwood (1982)[97]	1920s–1930s			0.04
Mulholland and Elwood (1982)[97]	1977–1979			0.3
Sarmiento and Sundquist (1992)[98]	1970s–1980s	0.8–0.9
Degens et al. (1991)[99]	1980s	0.7–0.8
Meybeck (1982)[100]	1970s–1980s	1.0
Dean and Gorham (1998)[101]	1970s–1990s			0.2
Meybeck (1993)[102]	1992	1.0		0.2
Aumont et al. (2001)[103]	1980s–1990s	0.8
Schlünz and Schneider (2000);[104] Aufdenkampe et al. (2011)[32]	1980s–1990s	0.8–0.9
Cole et al. (2007)[1]	1990s–2000s	0.9	0.8	0.2	1.9
Battin et al. (2009)[5]	1990s–2000s	0.9	1.2	0.6	2.7
Tranvik et al. (2009)[2]	1990s–2000s	0.9	1.4	0.6	2.9
Aufdenkampe et al. (2011)[32]	1990s–2000s		3.3
Regnier et al. (2013)[9]	1750	0.8	0.7	0.2	1.7
Regnier et al. (2013)[9]	2000–2010	1.0	1.2	0.6	2.7
Raymond et al. (2013)[8]	1990–2010		2.1
Deemer et al. (2016)[105]	1990–2010		2.7
Sawakuchi et al. (2017)[6]	1990–2016		2.5		5.1
Drake et al. (2018)[11]	1990–2016	1.0	2.5	0.6	5.1
this study	1900	0.9	1.6	0.4	3.0
this study	2000	0.9	1.9	0.5	3.3

The increase in C burial during the past century reported by existing studies (Table ) is well captured by DISC-CARBON (Table and Figure a). This is mainly the result of increasing C inputs (Figure a). Our estimated 0.5 Pg C yr–1 of C burial for the year 2000 is very close to 0.6 Pg C yr–1 for 1990–2016 reported in recent studies[2,5,9,11] and obviously higher than the earlier estimate of 0.2–0.3 Pg C yr–1 for the preindustrial era[9] and 1970s–2000s.[1,97,101,102] Our simulated C burial of 0.37 Pg C yr–1 in 1900 is higher than 0.04 Pg C yr–1 for the 1920s–1930s estimated by Mulholland and Elwood (1982)[97] because they only counted C burial in lakes and reservoirs and ignored the decline of the lake area since the 1920s–1930s when using the estimate for the 1970s as the basis to trace back the C burial during the 1920s–1930s. Our estimate of 1.9 Pg C yr–1 of the CO2 emission from global inland waters for the year 2000 is in the middle of the range of 0.8 Pg C yr–1 by Cole et al.[1] and 3.3 Pg C yr–1 by Aufdenkampe et al.[32] for the 1990s–2000s and is very close to the CO2-partial-pressure-data-based estimate of Raymond et al.[8] for 1990–2010 (2.1 Pg C yr–1). The estimated CO2 emission for the Amazon (Table ) is consistent with the estimate of 0.8 Pg C yr–1 by Rasera et al.[500], but twice the estimate of 0.5 Pg C yr.–1 by Richey et al.[600] The estimate of1.4 Pg C yr.–1 by Sawakuchi et al.[6] includes emissions from the many tidal floodplains in the lower Amazonian basin. This suggests that the lower estimates for global CO2 emissions shown in Table may be underestimates as they are close to the estimated emissions attributed to the Amazon River alone.

Our CO2 emission estimate for the year 2000 is somewhat lower than 2.5 Pg C yr–1 reported by Sawakuchi et al.[6] and Drake et al.[11] for the period 1990–2016. This temporal increase (between the periods 1990s–2000s and 1990–2016) and the increase from 0.7 to 1.2 Pg C yr–1 between the preindustrial era and the 2000s observed by Regnier et al.[9] are consistent with the increasing trend in the global freshwater CO2 emission from DISC-CARBON (Figure a). Our estimated 1.6 Pg C yr–1 of CO2 emission from inland waters for the year 1900 is higher than 0.7 Pg C yr–1 for the preindustrial era estimated by Regnier et al.[9] because of the increasing trend in the CO2 emission during 1750–1900 and the possible underestimation of historical emissions by Regnier et al.[9] due to ignoring the change in natural C inputs.

These increasing C inputs to the global inland waters have been accompanied by increasing C retention and in particular emission to the atmosphere, with the consequence that the river C export to the world’s oceans has remained stable (Figure a,b). Our simulated global TC export of 0.9 Pg C yr–1 in 2000 almost equals the estimates of 0.9–1.0 Pg C yr–1 for the period 1990–2016 in previous studies[1,2,5,9,11] (Table ). This lack of change in TC export was also observed in previous estimates of 0.8 Pg C yr–1 for the preindustrial era,[9] 0.7–1.0 Pg C yr–1 for the 1970s–1980s,[98−100] and 0.9–1.0 Pg C yr–1 for the most recent decade[11] (Table ).

Although our simulation results show that the global export of TC remained rather constant (in terms of quantity and forms), there are regional differences. For instance, C export into the Mediterranean Sea and the Black Sea showed a decline (Figure b). The river DOC export from Arctic rivers simulated by DISC-CARBON is 38 Tg C yr–1 for the period 1950–2000 (Figure b), which is in good agreement with the estimates of 25–36 Tg C yr–1 based on the observations in other studies.[94,96] Our simulated global DOC export of 129 Tg C yr–1 is ∼20% lower than the regression result by Global NEWS for the year 2000,[106] while our simulated global POC export exceeds the previous regression-based estimates by ∼50% for the 1990s.[28,30] One reason for the difference between our model and these regression-based estimates of POC export is related to the large uncertainties in the regressions used for calculating the total suspended sediments (TSS), which was the basis for their POC calculation. In this study, sediment dynamics is strongly linked with C cycling along the entire aquatic continuum and is thoroughly described in the DISC-CARBON module (Text SI1), which reduces the uncertainty.[52] Another reason is that the Beusen et al. estimate[28] is based on data for 19 European rivers, which may have induced a bias.

Overall, the temporal changes simulated by DISC-CARBON are consistent with the trends from earlier data-based and model-based assessments (Table ). In DISC-CARBON, the simulated long-term trends of the annual freshwater C fluxes indicate that global river basins have been balancing the increased inputs through more in-stream retention and emission to the atmosphere. The increasing C retention in inland waters may be closely related to the increasing global dam construction and thus increasing reservoir volume, which has increased from 7 to 3800 km3 parallel to dramatic land-use changes between 1950 and 2000 (Figure SI4). The result of the C budget from Regnier et al.[9] also shows similar increases in C input, emission, and burial and a relatively stable C export between the preindustrial era and the 2000s, and the increase in C burial is mainly attributed to that in reservoirs.

Future Improvements and Outlook

This study provides an integrated view and consistent quantification of the spatiotemporal changes in global freshwater C cycling over long time series, which can only be achieved with (i) a full and form-explicit process representation and (ii) spatially explicit and dynamic C inputs and environmental forcings. The long-term simulations of DISC-CARBON show good agreement with the available observations covering the second half of the 20th century (with close to 90% within a factor of 2, and no systematic errors), and the simulated trends for the full 20th century agree with the results of other studies based on various approaches.

However, our model assessment and sensitivity analysis demonstrate that improvement in some formulations and model input parameters could result in an improved description of the C fluxes in river basins and a better match of simulations with observations. Such an improved model would lead to not only a better understanding of the role of rivers in the global C cycle but also better information for policy makers about the influence of human interferences on emissions of greenhouse gases from inland waters to the atmosphere.

The sensitivity analysis clearly points to some parameters that have an overall influence on the model results. The first is the discharge, which affects several biological and physical processes important for C dynamics in river basins. To improve this, a better description of hydrology in low-order streams to replace the current parameterization will improve our C cycle model in headwaters. The HydroSheds dataset[107] is a good candidate to improve our PCR-GLOBWB model. The HydroLakes dataset[108] can be useful to improve the current data on lake and reservoir water volume.

Second, the sensitivity analysis pointed to the importance of C inputs from terrestrial systems and alkalinity and DIC inputs via groundwater. These C inputs are an important source of uncertainty in terms of their spatial distribution, organic/inorganic ratios, and forms (dissolved or particulate). For example, an improved estimate of terrestrial POC input from litterfall and its variation from headwaters to mainstreams is necessary to provide a more robust quantification and spatial estimate of CO2 emissions from freshwaters. Similarly, alkalinity input from groundwater is important but uncertain and can be improved by a better model for weathering and DOC input to aquifers and transformation to DIC. Furthermore, accounting for PIC dynamics will form a necessary extension of the model.[72]

Furthermore, a future major challenge to be tackled in global biogeochemical modeling frameworks is to include nutrient limitations to primary production (N, P, and silicon) and oxygen availability for respiration and C burial, particularly in reservoirs. Apart from advancing our capabilities to simulate freshwater C dynamics, this will allow extending the applications to other important issues such as global emissions of greenhouse gases (nitrous oxide and methane) from inland waters.

Finally, the current model simulates annual fluxes, while for many processes, it may be important to analyze fluxes at shorter time scales. A shorter time step (e.g., monthly) would allow for the simulation of seasonal C fluxes in river basins and export to coastal waters to better understand the impacts on coastal marine ecosystems. A first experiment for the Rhine River basin (Text SI3) showed that the model with annual and monthly settings yields similar results.