The importance of hydrology in routing terrestrial carbon to the atmosphere via global streams and rivers.

SignificanceStream/river carbon dioxide (CO2) emission has significant spatial and seasonal variations critical for understanding its macroecosystem controls and plumbing of the terrestrial carbon budget. We relied on direct fluvial CO2 partial pressure measurements and seasonally varying gas transfer velocity and river network surface area estimates to resolve reach-level seasonal variations of the flux at the global scale. The percentage of terrestrial primary production (GPP) shunted into rivers that ultimately contributes to CO2 evasion increases with discharge across regions, due to a stronger response in fluvial CO2 evasion to discharge than GPP. This highlights the importance of hydrology, in particular water throughput, in terrestrial-fluvial carbon transfers and the need to account for this effect in plumbing the terrestrial carbon budget.

The Earth’s water, carbon, and energy fluxes follow seasonal variations in the Earth’s solar radiation and climate variability (1, 2). As an integral part of terrestrial landscapes, streams and rivers receive significant water and carbon inputs from terrestrial and wetland ecosystems, which are further processed along the river to ocean continuum (3). As the largest carbon flux mediated by fluvial systems, carbon dioxide (CO2) emission from stream and river surfaces (4–7) is double the lateral carbon transport to oceans (8), yet its spatial and temporal variations are not fully resolved. Stream and river CO2 evasion changes considerably across space and time due to biogeochemical responses to climatic factors (3), the physics governing the transfer of gas across the water–air interface (9), and seasonal variations in the spatial extent of drainage networks (10, 11). However, seasonal variability of the flux has not been determined at the global scale, limiting our ability to understand controls at the macrosystem level.

The rate at which streams and rivers exchange CO2 with the atmosphere is determined by three factors: dissolved CO2 concentration (often expressed as an equivalent atmospheric partial pressure [pCO2]), water surface gas transfer velocity (k), and water surface area. To estimate flux at the monthly scale, all three factors need to be resolved at the same or finer temporal scale(s). To date, existing spatially explicit estimates of riverine CO2 emission at the global scale (4, 12) relied exclusively on pCO2 calculated from carbonate equilibria and historical archives of pH and alkalinity measurements. While these data have reasonable spatial coverage, the carbonate equilibria method is subject to inflated pCO2 estimates due to biased pH measurements (13) and alkalinity contribution from organic acids (14), particularly in low–ionic strength waters. These errors, although reducible within individual datasets (15), are difficult to correct for when scaling globally. This problem has significantly undermined calculations and understanding of the flux at the global scale. More importantly, although global estimates of the stream and river surface area and gas transfer velocity at mean annual discharge have been achieved (4, 16), their seasonal extent, a major driver of within-year variability of riverine CO2 flux, has not. This is largely because a temporally resolved reach-scale representation of global river hydrology has not been available until recently (17), and new understandings of aquatic surface area extent and water–air gas transfer rates are necessary to incorporate temporal variability into the riverine CO2 flux estimate.

We compiled a dataset of present-day direct pCO2 measurements in global streams and rivers from the literature. The dataset has 5,910 individual measurements of different months that cover all major freshwater ecoregions of the world (18), despite a small percentage (∼0.5%) of measurements from southern temperate rivers (). The dataset further has pCO2 measurements in all months from each freshwater ecoregion (open water months for the polar freshwater ecoregion) except oceanic islands and large river deltas that make up only 0.4% of the global land area (). These observations allowed for robust validation of the study’s results. Riverine pCO2 was statistically modeled against a set of watershed properties () in order to understand biogeochemical and geophysical controls on pCO2. Predictions of pCO2, k, and surface area were based on a new representation of the global river networks (the Global Reach-Level A Priori Discharge Estimates for Surface Water and Ocean Topography [GRADES] river networks) (17), which contains daily discharge estimates at ∼3 million individual river reaches over a 35-y period. Monthly CO2 flux estimates were achieved by coupling monthly pCO2 estimates driven by monthly watershed properties to monthly k and surface area estimates driven by the GRADES discharge. Spatial and temporal variability of the flux was finally investigated to demonstrate a strong modulation of the terrestrial (and wetland) carbon routing to the atmosphere via streams and rivers by hydrology.

Watershed pCO2 Controls

Understanding how stream pCO2 is related to different watershed properties is necessary for improving the mechanistic understanding of the flux and successful modeling of its variability across space and time. Below, we explore broad-scale watershed controls on river pCO2 via linear regressions with delineated watershed properties. Of the 26 watershed properties examined, soil respiration rates (autotrophic, heterotrophic, and total) explain the highest percentages of the pCO2 variability (coefficient of determination or R2 = 0.34 to 0.35) (). Soil organic carbon content (SOC), which measures organic carbon storage in soils, does not correlate with pCO2 (R2 = 0.01), in contrast to earlier studies that use SOC as a major geographical predictor for stream pCO2 (19). Gross primary productivity (GPP), which measures total carbon fixation by terrestrial plants as both biomass increase and plant maintenance respiration, shows a stronger correlation with pCO2 than net primary productivity (NPP), which measures terrestrial carbon fixation as biomass increase (R2 = 0.31 vs. 0.22). Temperature and precipitation are also correlated with pCO2 (R2 = 0.28 and 0.11, respectively), suggesting broad-scale climatic controls on stream pCO2. Wetland areas show a weak correlation with pCO2 at the global scale (R2 = 0.03) despite their importance in tropical lowland rivers (5, 20), suggesting probably local or regional influences in wetland-rich systems.

The fact that soil respiration rates are among the best predictors of river pCO2 highlights close linkages between soil carbon dynamics (in particular, soil CO2 stripping) and watershed carbon loss through water surface CO2 evasion, widely recognized in headwater streams (21). Watershed properties suggestive of ecosystem carbon storage (e.g., SOC) or net photosynthetic fixation (e.g., NPP) are, however, comparatively less relevant. Considering that direct terrestrial–stream linkages are more significant in small headwater systems (21) or at terrestrial–riverine interfaces in larger systems (6), the strong predictability of stream pCO2 by soil respiration or terrestrial GPP (R2 = 0.31) is also likely caused by broad-scale geographical and seasonal synchronicities between these fluxes (2, 22) and an ensemble of riverine CO2-relevant hydrologic, biogeochemical, and biospheric processes (3) driven by the same climatic factors. Most importantly, temperature and precipitation are two common factors driving terrestrial carbon fluxes and exchanges (1, 2), which underlie the supply of reduced or respired contemporary terrestrial carbon to fluvial systems. Climate factors also strengthen instream processes by enhancing organic carbon metabolism under higher temperatures and/or strengthened hydrologic connectivity (3).

Watershed slope and elevation are negatively correlated with river pCO2 (R2 = 0.33 and 0.19, respectively). Although slope and elevation are partially correlated with each other (e.g., steep terrain is found more often in higher elevations), they point to distinct watershed controls on the pCO2 variability. While high elevations feature colder climates, low terrestrial productivity, and soil respiration rates, the negative effect of watershed slope is also related to rapid release of CO2 from high-gradient, turbulent water columns (19, 23). Considering that pCO2 in flowing waters reflects an equilibrium between the rate of source inputs and surface evasion, the effect of watershed slope is significant as it manifests the physical and geomorphological controls of river pCO2 via surface evasion (9).

pCO2 Modeling and Seasonal Variability

For a dataset of mixed spatial and temporal observations, it is essential that a model predicts accurate pCO2 across both space and time. Below, we examine how modeling stream pCO2 against watershed properties resolved at the monthly scale vs. modeling them against annual watershed properties differs in predicting the spatial and seasonal variabilities of pCO2. The modeling was done with a random forest (RF) regression model, which accounts for nonlinearities and shows a much stronger predictability than a reduced multilinear (RML) regression model (R2 = 0.77 vs. 0.38) (). We show that modeling pCO2 against monthly watershed properties yields much stronger predictability for pCO2 than modeling pCO2 against annual values of watershed properties in each individual month (R2 = 0.6 to 0.87 vs. R2 = 0.29 to 0.58) (Fig. 1 and ). Furthermore, comparing predicted pCO2 against site-level direct seasonal observations suggests that the monthly watershed properties vs. pCO2 approach yield much more reliable seasonal variations than the annual watershed properties vs. the pCO2 approach (Fig. 1 and ). This analysis suggests that unavoidable cross-month averaging of watershed predictors and pCO2 by the annual model results in loss of essential seasonal information critical for seasonal variability prediction.

Fig. 1.

Modeling stream pCO2 against monthly watershed predictors using an RF model. (A) Modeled and measured pCO2 values show high correlations (R2 = 0.6 to 0.87) in separate months. (B) Predicted pCO2 (normalized to site averages) shows seasonal variations agreeing well with direct observations across climate regions. (C) Histograms of the direct pCO2 measurements by climatic regions. Numbers indicate the number of measurements in each climatic zone. Dashed vertical lines indicate the atmospheric pCO2 (∼380 µatm).

Our estimate shows high consistency with field pCO2 measurements worldwide (). On the whole, pCO2 ranges from below the atmospheric average (380 µatm) to over 10,000 µatm at the reach level (Fig. 2). Spatially, high pCO2 (e.g., >3,000 µatm) is found in hot/humid lowlands, including the tropical Amazon [e.g., reported pCO2 is 3,000 to >12,000 µatm (6)], central Congo [e.g., reported pCO2 is 300 to 17,000 µatm (5)], Southeast Asia, southern subtropical United States, and India. Low pCO2 (e.g., <1,000 µatm) is found in frigid/arid climates where terrestrial inputs and in situ production are likely much weaker (24). Low pCO2 is also found in elevated, steep-terrain areas where stream CO2 losses due to surface evasion are much stronger (19). As an example, low pCO2 (5th to 95th quantiles: 590 to 1,310 µatm; close to the median of 860 µatm reported in ref. 25) of the upland Tibetan Plateau rivers could be a result of both weak inputs and strong evasion in this area. Across the climate zones, average pCO2 in tropical rivers is 40 to 70% higher than in temperate and Arctic rivers (i.e., 2,560 vs. 1,540 to 1,810 µatm) (Fig. 2 and ), which is consistent with the broad-scale geographical pCO2 variability reported by major regional studies (5–7, 26).

Fig. 2.

Maps showing spatial distribution (Left) and monthly variations (expressed as C.V. of monthly values; Right) of pCO2 (A and B), gas transfer velocity (k; C and D), surface CO2 efflux (E and F), and surface area extent (G and H) of global streams and rivers. C.V. was calculated from surface area–weighted mean monthly values in each climatic zone (). Dashed lines (B, D, F, and H) indicate latitudes that separate climatic zones. Note that river surface area is shown as a percentage of land area per HydroBASINS Level 04 basin (53).

Monthly pCO2 variability is greater in Arctic and temperate regions than in tropical regions (coefficient of variation [C.V.]: 8 to 11 vs. 5%) (Fig. 2), in alignment with stronger climatic variability in these areas (i.e., much larger seasonal temperature and precipitation differences) (27). In temperate and Arctic rivers, the highest pCO2 is generally found in the midsummer of each hemisphere (i.e., July and February in the Northern and Southern Hemispheres, respectively), which is 1.2 to 1.4 times the lowest pCO2 in midwinter (Fig. 1). Although in situ plant drawdown is greater during summertime, this is not large enough to drive monthly patterns (28). This seasonal difference is within the range (one to three times) often reported in temperate streams and rivers (28–30). Tropical regions show higher pCO2 in flooded seasons (Fig. 1), most likely driven by variable hydrologic connectivity to wetlands, a significant CO2 source to tropical rivers (20).

Gas Transfer Velocity

Monthly gas transfer velocity was estimated from reach-level slope (S) and flow velocity (V), which together correspond to the decaying energy dissipation along river networks (k600 = 2,841SV + 2.02) (9). We used a gauge-derived discharge to velocity (Q-V) relationship (i.e., the US Geological Survey equation in ref. 4) to estimate global flow velocities. Comparing predicted flow velocities against those derived from hydrologically routed flows from local and national river networks () suggests that the low-flow biased Q-V relationship derived from the stream gas tracer experiment (9) significantly underestimates flow velocity, particularly for small rivers, which are significant to the global river network surface area (e.g., 0.26 vs. 0.1 m s−1 in rivers of <30 m3 s−1; these rivers make up >90% of the 3 million GRADES river reaches). The gauge-derived Q-V relationship, however, predicts reliable flow velocity over a broad range of discharge (0.01 to 20,000 m3 s−1). In comparison with flow velocity, slope measurements by various digital elevation models (∼30-m resolution) are relatively reliable () and not expected to cause high uncertainty. Predicted k is in good alignment with those reported from regional studies ().

Bubble-mediated gas exchange has been recently recognized as an important mechanism in high-energy streams (23). This nonlinear mechanism is caused by strong bed friction and significant in mountain streams of shallow water depth and high bed roughness, for which steep terrains act as an integral driver (23). Using a 0.01 (unitless) “mountainous” slope cutoff for streams that have significant bubble-mediated gas exchanges (), high-slope mountain streams make up 14% of the total GRADES river reaches. For these high-energy systems and affected extrapolated areas, k was estimated according to a suggested power law relationship with stream dissipation energy (23). Predicted k for mountain streams is close to that reported in Horgby et al. (19) (16 to 23 vs. 26 m d−1) ().

The highest k (e.g., >16 m d−1) is found in major elevated regions of the Himalayas, the Rockies, the Andes, the Alps, and eastern Africa, and the lowest k (e.g., <3 m d−1) is found in major lowlands of low to moderate humidity, including the great plains of North American and western Siberia (Fig. 2). In alignment with decaying energy dissipation rate (eD) along stream order (9), k decays significantly with system size (). Monthly variability in k is much stronger in Arctic and northern temperate rivers than in tropical and southern temperate rivers (with C.V. of 21 to 43 vs. 2 to 4%) (Fig. 2) because of the strong annual hydrological variability in northern regions of strong continental climates (27, 31). The highest gas transfer rates are found in midsummer of each hemisphere (), largely because of the higher late spring/summer flows widely found in many climate types [e.g., flow regime six to eight (31) in monsoon- and ice melt–affected climates] ().

CO2 Efflux and Seasonal Variability

Predicted CO2 efflux is in general agreement with that reported from regional studies (). Tropical rivers have much higher CO2 efflux than Arctic and temperate rivers (3,220 vs. 1,750 to 2,280 g C m−2 y−1) (Figs. 2 and 3). CO2 efflux in temperate rivers is also close to that reported by Butman and Raymond (7) for the northern temperate region (2,060 vs. 2,370 g C m−2 y−1). Seasonally, similar to pCO2 and k, Arctic and northern temperate rivers exhibit much greater monthly variability in CO2 efflux than tropical and southern temperate rivers (C.V.: 25 to 37 vs. 4%, respectively) (Fig. 2). For instance, the highest monthly average CO2 efflux is twice the lowest average in northern temperate rivers (2,750 vs. 1,360 g C m−2 y−1), highlighting strong seasonal variability in CO2 flux from these rivers.

Fig. 3.

Monthly surface CO2 efflux in global streams and rivers. Ice-covered regions are grayed out for each month.

Monthly River Surface Area

Monthly river surface areas were estimated based on a combination of downstream hydraulic geometry (DHG) and at-a-station hydraulic geometry (AHG) and an improved land surface area segregation scheme based on watershed hydrology ( and ). Total annual surface area of global streams and rivers is 811,000 km2, close to the remote sensing–based estimate in Allen and Pavelsky (16) (105%, 773,000 km2) but 30% higher than the raw surface area in Raymond et al. (4) (624,000 km2). Ephemeral area or surface area loss due to temporary drying up of intermittent rivers, estimated using two different methods (), is 4 to 7% of the total river surface area at the global scale and related to watershed runoff. For example, low watershed runoff is responsible for the high ephemeral extent (13 to 16%) found in southern temperate rivers where arid Australian central and western rivers (with watershed runoffs of <50 mm y−1) make up a significant portion of the total surface area. Surface area loss due to winter ice coverage is much more variable across climate zones and months than surface area loss due to drying, despite a rudimentary accounting for ice dynamics in this analysis (Fig. 3 and ). Ice coverage is most prevalent between October and April, ranging from 53 to 92% and from 5 to 40% in Arctic and northern temperate basins, respectively. Globally, ice coverage extent ranges from 0 to 27% in different months, lower than the 0 to 56% ice coverage from a recent estimate (32), probably because a rather conservative temperature cutoff (−4 °C) was used for ice coverage initiation. Overall, ice coverage causes a 12% loss to the global river surface area. Accounting for drying up and ice coverage, the monthly area of flowing waters ranges from 500,000 to 854,000 km2 and averages 672,000 km2 (), which is 25% higher than the estimate in Raymond et al. (4). The analysis suggests strong cross-month changes in river surface area caused by seasonal variability in global hydrology and ice coverage (Fig. 3). In particular, monthly river surface area is different by a factor of >2.3 in northern temperate watersheds (C.V.: 30%), and that in Arctic regions drops to ∼0 in frozen seasons (C.V.: 90%) (Fig. 2), showing much stronger variability than tropical and southern temperate watersheds (C.V.: 2 to 9%) (Fig. 2).

CO2 Emission and Seasonal Variations

Annual CO2 emission from global streams and rivers amounts to 2.0 ± 0.2 Pg C y−1, including corrections for an enhanced release of ∼60 Tg C y−1 from spring ice melting in northern rivers ( and ) and a continued release of ∼50 Tg C y−1 from dried river channels of intermittent rivers (33). Our estimate is in general agreement with regional estimates from the literature (). An exception is wetland-dominated tropical lowland systems (i.e., the Congo and the Amazon), where our estimate is >20% lower. We suggest that our estimates are conservative because the methods were only able to account for emissions from flowing channels, and those from fringing floodplains or wetlands were not included (6). This becomes evident when comparing our estimates with separate CO2 emission estimates for the Amazonian streams or rivers (floodplains and wetlands excluded) (34), where our estimates are comparable or larger (). Higher discrepancy between the whole-basin estimates, however, warrants clearly defined boundaries and fully resolved temporal estimates devoted specifically to these systems (20).

Our estimate is slightly higher (∼11%) than the earlier estimate of 1.8 Pg C y−1, suggesting that the effect of reduced pCO2 from resorting to direct measurements was largely compensated for by higher gas transfer rates and water surface area, particularly from the smallest streams. Across climate regions, tropical rivers are responsible for 57% of the global emission, more than temperate and Arctic regions combined (30 and 13%, respectively), suggesting a dominant role of tropical rivers in the flux. The climatic zone estimates are in good alignment with earlier estimates from geographical extrapolation of regional surveys. For instance, total riverine CO2 emission for northern temperate rivers was estimated to be 0.54 Pg C y−1 (7), in comparison with 0.49 Pg C y−1 in this analysis, further corroborating the robustness of the analysis.

Seasonally, CO2 flux varies from 112 to 209 Tg C mo−1 (or 1.35 to 2.54 Pg C on an annual basis) (), suggesting a difference in monthly emissions of a factor of approximately two. The highest emission is found from May to August (203 to 209 Tg C mo−1), and the lowest is from November to February (112 to 115 Tg C mo−1), corresponding to the seasonality of the Northern Hemisphere where 75% of the land surface is located. Monthly variability is most prominent in Arctic and northern temperate rivers. In northern temperate rivers, the emission varies across months by a factor of 4.5, and in Arctic rivers, the emission drops to ∼0 in frozen seasons and rises until the end of spring freshets (i.e., June) (Fig. 3 and ). The large seasonal variations suggest significant temporal changes in the flux and unreliability of using a single annual estimate (4–7) for the seasonally varying fluvial fluxes (35). While seasonal pCO2 changes are ubiquitous across climate regions (C.V.: 5 to 11%) (Fig. 1), changes in k and water surface area are progressively more significant toward northern colder regions (Fig. 2 ), for which seasonal variations in watershed hydrology are a key factor (31). It is emphasized that seasonal watershed biogeochemistry and physical constraints imposed by seasonal changes in water surface area and surface water turbulence are important dynamics for seasonal riverine greenhouse gas evasions.

Terrestrial Carbon Routing to the Atmosphere Modulated by Water Throughput

Stream and river CO2 emissions are driven by terrestrial carbon inputs, either as the stripping and delivery of CO2 from soils (21) or flooded lands (20) or the dissolution/erosion of terrestrial or wetland organic carbon that is later oxidized within drainage networks (3) (). Recent work has argued that stream and river CO2 evasion is ultimately balanced by terrestrial GPP (36). Predominant young ages of fluvially evaded CO2 support a strong coupling between the fluvial emission and contemporary terrestrial GPP rather than old terrestrial carbon stocks (37, 38), although exceptions may occur (39). Globally, fluvial CO2 emission makes up ∼1.8% of global GPP (109 Pg C y−1) (22) despite a lower river to land surface area ratio (0.5%), reinforcing the concept of streams and rivers as hot spots for terrestrial–atmospheric carbon exchange. Here, we demonstrate that this emission percentage is highly variable, ranging from negligible (<0.2%) to as high as >3.6 to 18% across regions (Fig. 4). Higher percentages are found in both the humid tropics, where the riverine emissions are high, and high latitudes, where terrestrial productivity rates are low and surface soil organic carbon stocks are high (40).

Fig. 4.

Terrestrial carbon routing to the atmosphere modulated by water throughput. (A) Map showing stream CO2 emission as a percentage of terrestrial GPP. (B) Stream CO2 emission as a percentage of GPP increases in watersheds of higher water throughput (i.e., discharge). (C) The emission percentage increases linearly to logarithmic discharge for watersheds with discharge of >100 mm y−1 and annual air temperature of >8 °C. (D) Land area normalized stream CO2 emission rate (CO2e) and (E) the emission ratio (expressed as a percentage of GPP) scale closely with the river to land area ratio (expressed as percentage of river surface area). F and G show the power law relationships fitted between GPP and watershed discharge (F) and between the emission rate (CO2e) and watershed discharge (G), respectively. The relationships suggest a steeper increase to discharge for the emission rate than for GPP (slopes: 0.67 vs. 0.39). In B, bars are color coded by watershed discharge levels; error bars indicate SD within each watershed group. Data points in C, F, and G are color coded by mean annual air temperature (degrees Celsius). The map and all relationships are based on HydroBASINS Level 04 (53).

Part of the variability in the percentage of terrestrial GPP routed to the atmosphere via streams and rivers can be explained by water throughput or discharge (Fig. 4). Terrestrial fluxes are rarely looked at through the lens of river discharge, despite discharge being one of the major terms in terrestrial hydrology (41). Here, we find that after we get above a discharge (normalized to watershed area, same below) of ∼100 mm y−1 (which occurs on ∼70% of the global watershed surface), there is a clear relationship between discharge and the ratio of stream CO2 evasion to GPP across watersheds (Fig. 4). The relationship is impacted by a number of high-latitude watersheds flanking the Arctic Ocean (Fig. 4) and becomes more linear when the analysis is limited to watersheds with mean annual temperatures of >8 °C (emission percentage = 0.27lnQ + 0.11, R2 = 0.24) (Fig. 4). The relationships between precipitation and evapotranspiration and this ratio are much weaker (R2 = 0.07 to 0.11) ().

We posit two major mechanisms leading to a greater percentage of terrestrial GPP being shunted to streams and rivers for atmospheric evasion at higher discharge. The first is an increase in the areal extent of drainage networks at higher discharge, demonstrated by a close cross-region correlation between the percentage of river surface area, emission rates, and the emission ratios (Fig. 4 ). Within-watershed studies have also found a higher river surface area during times of greater discharge (10, 42). Stream and river scientists argue that increases in the degree of connectivity between terrestrial and fluvial systems are an important consideration in the transfer of terrestrial constituents to drainage networks (10). We argue that greater connectivity in wetter regions and during wetter months leads to a higher amount and fraction of terrestrial GPP being laterally transported (Fig. 4 ).

The second mechanism is the differing responses of fluvial vs. terrestrial fluxes to water throughput. Water is a fundamental constraint on terrestrial plant carbon fluxes. In brief, in order to achieve high rates of CO2 uptake through the stomate (i.e., GPP), terrestrial plants must transpire precipitation water before it enters the river network (43). Plants, however, have a finite need for water to balance GPP due to other limitations on GPP, such as nutrients or light (44, 45). So, we might expect that watersheds with high precipitation have “excess precipitation” that is either stored locally or transported to become river discharge (41, 46). Although quantifying vegetation response to water availability at large spatial scales is difficult, terrestrial studies have demonstrated a response of terrestrial carbon fluxes, including GPP, to the magnitude of precipitation (44, 47). Here, we report a nonlinear relationship between terrestrial GPP and discharge (lnGPP = 0.39lnQ + 4.7, R2 = 0.66) (Fig. 4).

The nonlinear relationship suggests a saturation of GPP per increase in discharge (i.e., lower ΔGPP per ΔQ) at high precipitation, in line with terrestrial studies that demonstrate a lower rain use efficiency for terrestrial plants in wet climates (44, 47). The increase in stream CO2 evasion with discharge, however, does not show a strong saturation and has a steeper response to discharge (lnCO2 emission rate = 0.67lnQ − 1.2, R2 = 0.84) (Fig. 4). This is in line with stream studies that demonstrate an increase in CO2 evasion and terrestrial carbon export at higher discharge (48, 49). Thus, stream CO2 evasion increases faster than GPP with increasing discharge (slope: 0.67 vs. 0.39), leading to a change in the ratio of fluvial CO2 evasion to terrestrial GPP, with a greater percentage of GPP being laterally delivered to and evaded from drainage networks in wetter watersheds (Fig. 4).

A larger ratio of stream/river CO2 flux to terrestrial GPP in regions of higher discharge provides initial evidence that the global terrestrial water cycle determines in part the relative importance of connectivity between terrestrial and fluvial ecosystems. Recent work has argued that integrating lateral fluxes into terrestrial carbon budgets is important to the estimate of soil respiration (35). These results indicate that improvement can be made to the plumbing of the terrestrial carbon budget if models can account for a differential response of terrestrial–atmosphere vs. terrestrial–drainage network carbon fluxes to water throughput. Furthermore, cold climates with moderate to high watershed discharge and large soil organic carbon stocks might be transporting as much as 5 to 18% of terrestrial GPP (Fig. 4 ) to drainage networks in the form of soil CO2 and organic carbon, highlighting the need for understanding carbon dynamics in these areas.

Materials and Methods

The materials and methods are described in complete detail in . The central features are summarized as follows.

To estimate monthly CO2 emissions from global streams and rivers, the three factors (pCO2, k, and surface area) determining aquatic CO2 fluxes need to be resolved at the same (i.e., monthly) or finer temporal scale(s). GRADES river networks (17) are a vector-based global drainage network dataset derived from the 90-m Multiple Error-Removed Improved-Terrain Digital Elevation Model (50). GRADES contains ∼3 million river reaches and continuous daily discharge estimates at these river reaches over a 35-y period (1979 to 2014). In this analysis, the GRADES river networks were used as the underlying hydrographic infrastructure for global river network CO2 emission estimates. Monthly pCO2, k, and surface areas (raw, ephemeral, and ice covered) were estimated at each GRADES river reach. Estimated river pCO2 values were driven by monthly watershed predictors (e.g., soil respiration, GPP, etc.), while k and surface areas were driven by monthly discharge from GRADES.

A dataset containing 5,910 direct fluvial pCO2 measurements was compiled from the literature (Dataset S1). We excluded pCO2 calculated from pH and alkalinity or measurements from open waters, considering that these values are subject to calculation errors (15) or affected by different dynamics than in streams and rivers (30). To predict pCO2 at GRADES river reaches, the pCO2 values were modeled against a set of watershed properties (climatic, geomorphic, and terrestrial carbon-cycling related) using an RF model and an RML regression model (R program, v3.6.2). The RF model outperformed the RML model (as suggested by lower model residuals and less significant bias along the pCO2 range) () and was selected for reach-level pCO2 prediction. Furthermore, we demonstrate that modeling individual pCO2 values against monthly watershed properties yields much more reliable seasonal variations than modeling site-average pCO2 against annual watershed properties (). The RF model was thus established using individual pCO2 values and monthly watershed predictors. Monthly pCO2 values at GRADES river reaches were then estimated using monthly watershed predictors at the river reaches as inputs.

GRADES river reach gas transfer velocities (k) were estimated from channel slope (S) and flow velocity (V; k600 = 2,841SV + 2.02, k600 is the gas transfer velocity normalized to a Schmidt number of 600), which together correspond to the decaying eD along river networks (9). GRADES provides channel slope estimates, which are in high agreement with other products (e.g., National Hydrography Dataset Plus) (). To scale flow velocity, we coupled monthly discharge from GRADES to a gauge-derived Q-V relationship (i.e., lnV = 0.12lnQ – 1.06) (4), which predicts reliable V over a broad range of river discharge (0.01 to 20,000 m3 s−1) (). For streams affected by high bubble-mediated gas exchanges in steep terrains, k was estimated using a reported power law relationship between k600 and the eD (lnk600 = 1.18lneD + 6.43) (23). A slope cutoff of 0.01 was used to differentiate low- vs. high-energy alpine streams ().

Monthly river reach widths were estimated by combining DHG and AHG relationships for width (). The DHG relationships were established at the mean annual discharge, which coupled to mean annual discharge from GRADES, were used to scale river widths along the GRADES river networks. AHG relationships, coupled to monthly discharge from GRADES, then allowed for scaling temporal variations in river widths based on widths estimated at the mean annual discharge. AHG width exponents at GRADES river reaches were predicted using a multilinear relationship established between flow characteristics (e.g., C.V. of daily discharge and runoff) and the exponent (). Reach length was provided by the GRADES river networks. Surface areas were calculated from monthly widths and reach length.

The GRADES river networks start channelization at ∼25 km2 (with first stream-order widths of ∼2 to 6 m) and miss the smallest streams and rivers. For these streams, k and surface areas were estimated by extrapolating related river hydromorphological or hydraulic characteristics (i.e., width, length, and k) according to stream-order scaling laws (51). We used a width of 0.3 m as the cutoff for extrapolation, which is the mean of median headwater stream widths found in a field study (52). pCO2 at the extrapolated stream orders was assumed to be the same as that in first-order streams of the GRADES.

It is essential to group land areas with similar flow characteristics together for extrapolation considering that the number of extrapolated stream orders varies spatially and temporally depending on watershed hydrology. An improved hydrology-based segregation scheme was, therefore, employed, which separates the global land surface area into a total of 78 cross-region basins based on watershed runoff levels and flow regimes (i.e., monthly flow distribution) (31) (). The segregation was also based on the HydroBASINS Level 04 watersheds (53). The cross-region basins were used as the basic land surface unit for drainage network surface area and CO2 emission estimates.

Ephemeral surface areas were corrected for using a nonlinear relationship developed between the mean annual discharge and the ratio of days a stream/river is dry (i.e., zero discharge) within a year (). Ephemeral surface areas for extrapolated stream orders were corrected for by scaling and extrapolating the estimated ephemeral extents of the GRADES river networks by stream order. Ice-covered river surface areas were corrected for by identifying land surfaces that are under −4 °C for each month (Fig. 3). CO2 emissions during ice-melt periods were accounted for by applying a reported ice-melt emission ratio (i.e., 0.17) (54) to flowing-water CO2 emissions estimated for each latitudinal band and month ().

Uncertainties associated with the monthly CO2 emission estimates were evaluated by considering two major sources of error in the current analysis: pCO2 error from the RF model and discharge error from GRADES (). A Monte Carlo process was pushed through all steps of the CO2 emission estimation procedures, and uncertainties were calculated as the 1δ deviation of the simulated emission magnitude distributions.