Hydrologic controls on pCO2 and CO2 efflux in US streams and rivers

Streams and rivers emit petagrams of CO2, yet there is little known about how discharge variability impacts stream CO2 at broad scales. This article analyzed changes of pCO2 and CO2 efflux with discharge in 813 conterminous US streams and rivers. Half of the streams and rivers (49%) showed positive pCO2–discharge responses, while significantly more (80%) showed positive CO2 efflux–discharge responses. More CO2 efflux therefore occurred at higher discharge. While CO2 efflux tended to increase with discharge at all stream orders, pCO2 tended to decrease with discharge in small rivers (stream order 1–6) and increase in large rivers (stream order 7–10). Small streams and rivers tended to release more CO2 at higher discharge due to highly responsive gas transfer velocity. While high CO2 efflux in small streams (stream order 1–3) pointed to strong terrestrial input, high CO2 responsiveness in large rivers was likely due to increased riverine metabolism and river–corridor exchanges.

Discharge (Q) is a core attribute of river networks. Concentration-discharge analysis contains information on production, spatial distribution, and watershed-hydrology interactions and is an important tool for mechanistic analysis of transport of suspended particulates, weathering solutes, nutrients, and dissolved organic matter by river networks (Lefrançois et al. 2007;Moatar et al. 2017;Diamond and Cohen 2018). Concentration-discharge relationships also form the base for flux-discharge integrations for estimating riverine fluxes (Raymond and Saiers 2010;Wilson et al. 2013) and therefore are important in the context of climate change where shift in precipitation regimes is a core component.
Streams and rivers are globally significant CO 2 emitters to the atmosphere Borges et al. 2015), where both the concentration (measured often in equivalent partial pressure of dissolved CO 2 or pCO 2 ) and surface efflux of CO 2 are thought to be strongly affected by Q. Although *Correspondence: shaoda.liu@yale.edu Author Contribution Statement: PAR conceived, designed, and guided the analysis. SL conducted data analysis, constructed tables and illustrations, and wrote the manuscript. Both authors contributed to interpretation of results and revising the draft.
Date Availability Statement: Water-quality measurements, stream order and channel slope information, and hydraulic equations of conterminous USGS sites used for this analysis can be found in the Environmental Data Initiative repository at https://doi.org/10.6073/pasta/ 5ed5b7b2a82f2381e759d519f27ff865.
Additional Supporting Information may be found in the online version of this article.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
there are scattered field studies suggesting increase or decrease of either pCO 2 (Dinsmore et al. 2013;Almeida et al. 2017) or CO 2 efflux (Long et al. 2015;Liu et al. 2017) with Q, a systematic regional to continental scale analysis is not available. Understanding the role of discharge however is essential for both validating the scalability of field observations and revealing relevant environmental controls on hydrologic responses of CO 2 concentration and emission. Furthermore, how changing pCO 2 and gas transfer velocity (k) interplay to control CO 2 efflux is a topic rarely researched in previous studies (Long et al. 2015;Liu et al. 2017).
Unlike nongaseous solutes, the transport of CO 2 by rivers is separated into two parts: downstream transport (the pCO 2 or concentration term) and surface evasion (the efflux term), the latter making up a significant carbon flux in rivers of all sizes and spatial scales (Butman and Raymond 2011;Raymond et al. 2013;Wallin et al. 2013;Borges et al. 2015). The two parts are inherently dependent on each other, for instance, while concentration promotes evasion, evasion depletes CO 2 . Although there is a mutual dependency, the two processes do have different drivers. While the concentration is determined by rates of source replenishment and evasion, efflux depends on both concentration and k, the latter being a function of surface-water turbulence and river hydraulics (Alin et al. 2011;Raymond et al. 2012). Understanding how pCO 2 and CO 2 efflux respond differently to Q therefore represents a new and important topic with respect to improving regional and global models of river CO 2 evasion.
We used a water quality dataset from the US Geological Survey (USGS) to quantify pCO 2 and CO 2 efflux in US streams and rivers. We researched the relationships between pCO 2 , CO 2 efflux, and Q and how the hydrologic responses of stream CO 2 were related to Strahler stream order of the river network and connectivity to terrestrial systems. This is the first systematic quantification of CO 2 -Q responses in streams and rivers across a broad spatial extent to our knowledge.
We coupled the water quality dataset with National Hydrography Dataset Plus (NHDplus) (http://www.horizon-systems. com/nhdplus/), where stream segment slope and stream order information are available for central flow lines of all streams and rivers in the United States. We further coupled the dataset with a watershed attribute dataset previously delineated by Raymond (2017), where drainage areas were available and used to normalize Q to water yields per unit area. A total of 494,549 measurements from 4818 sites were left with paired slope, stream order, and drainage area information.
A separate dataset containing coupled hydro-geometric measurements of stream/river Q, width (w), velocity (v), depth (h), and cross-sectional area was available from the USGS NWIS's surface water web interface (https://waterdata.usgs. gov/nwis/sw). This dataset was used to derive separate hydraulic geometry relationships (equations below) for each site, which was also paired with the water-quality measurements.
where a, b, and c are the hydraulic coefficients and d, e, and f are the hydraulic exponents.
pCO 2 overestimation and data binning One problem of calculating pCO 2 from raw pH and alkalinity measurements is that the method yields overestimated pCO 2 values at low pH, hampering correct interpretation of CO 2 trends. These overestimations result from dysfunctioning of the carbonate equilibria in acidic waters of high-organic alkalinity (Abril et al. 2015) or simply erroneous pH measurements especially in waters of low buffering capacity. In order to avoid these pCO 2 overestimations, we employed a binning process at each site. According to the binning method, water-quality measurements of the same site were arranged in a descending order according to Q. The measurements were then segregated into bins of size 10. Within each data bin, the median value was calculated and used for pCO 2 and CO 2 efflux calculation and hydrologic responses analysis (see the method sections that follow). We constrained total number of measurements to be ≥ 180 at each site. This ensured a minimum of 18 points at each site at a binning size of 10 and yielded the final dataset of 303,583 measurements from 813 sites for the current analysis ( Fig. 1). A summary of geophysical and hydrological properties of the stream/river sites and pCO 2 according to stream order is available (Supporting Information Table S1).
A sensitivity test (Supporting Information Figs. S1, S2) was performed in order to examine the effect of binning size on effectiveness of screening out extremes and preserving trends. It was found that a binning size of 10 was able to avoid most pCO 2 overestimations (as indicated by quick decrease of maximum pCO 2 values with binning size, Supporting Information Fig. S1) and preserve regression trends at the same time (Supporting Information Fig. S2). A binning size of 10 was therefore adopted in this analysis.

Liu and Raymond
Hydrologic controls on river CO 2 efflux Evasion flux calculation Water pCO 2 was calculated according to freshwater carbonic acid dissociation equilibria along with dissociation constants from Clark (1997). CO 2 efflux (in g C m −2 yr −1 ) was estimated as the product of dissolved CO 2 concentration gradient at the water-air interface and gas transfer velocity (k): where [CO 2 ] w and [CO 2 ] a represent CO 2 concentration (mmol L −1 ) in water and in equilibrium with a constant atmospheric pCO 2 of 390 μatm, respectively, and k is in m d −1 . A temperature-sensitive Henry's Law constant (H, in atm × L/mol) was used to convert partial pressures to dissolved CO 2 concentrations: − log 10 H ð Þ ¼ − 7 × 10 5 T 2 + 0:015T + 1:11 ð6Þ k normalized to a Schmidt number of 600 (k 600 ) was estimated using stream velocity (m s −1 ) and slope (S, unitless) (Raymond et al. 2012): where S was adopted from the NHDplus dataset and v calculated from site-specific hydraulic geometry equations and gauge Q (Eq. 3). k was converted from k 600 using Schmidt number at median temperatures of each site with the following equation (Raymond et al. 2012): The Schmidt number for CO 2 at different temperatures ( C) in freshwater was calculated following where 600 is the Schmidt number of CO 2 in freshwater at 20 C (Raymond et al. 2012).

Hydrologic response modeling
In order to model changes of variables (pCO 2 , k, and CO 2 efflux) with Q, we transformed median measurements and derivative variables to their natural logarithmic values. This relieved concerns on non-normality of the raw data but also allowed description of the CO 2 -Q dependencies with power law relationships similar to those of river hydraulic geometry (Eqs. 1-3): pCO 2 ¼ e ic Q Sc , k ¼ e i k Q S k , and CO 2 efflux ¼ e iF Q SF . Or, where S c , S k , and S F are linear regression slopes between ln pCO 2 , ln k, ln CO 2 efflux, and ln Q, respectively, and i c , i k , and i F the y-axis intercepts. The regression slopes (S c , S k , and S F ) were used to indicate directions (positive or negative) and magnitude of responses of relevant variables to Q. This model was not able to model changes with Q when CO 2 efflux was negative (which makes up 1.5% of the total calculated effluxes) due to the logarithmic transformation.

Results and discussion
About half (49%) of the sites were associated with positive log-linear slopes or increasing pCO 2 with Q (Fig. 2). Gas

Liu and Raymond
Hydrologic controls on river CO 2 efflux transfer velocity (k) responded positively to Q at all sites (with log-linear slopes of 0-0.83), in line with physical controls on gas transfer velocity by stream hydraulics and surface-water turbulence (Zappa et al. 2007;Alin et al. 2011;Raymond et al. 2012), which both often scale positively with Q. Thus, CO 2 efflux increased with Q in a significantly higher proportion (80%) of the US streams and rivers, suggesting that the majority of streams and rivers have higher CO 2 efflux at higher Q. Three main types of streams and rivers were identified ( Fig. 1): (1) streams/rivers with both positive concentration and CO 2 efflux responses, (2) streams/rivers with negative concentration but positive CO 2 efflux responses, and (3) streams/rivers with both negative concentration and CO 2 efflux responses, making up 48%, 32%, and 20% of the total stream and river sites, respectively (Supporting Information  Table S2). Among sites with positive pCO 2 responses, the majority (99%) also showed positive CO 2 efflux responses, suggesting a positive concentration response largely determined an increasing CO 2 efflux-Q response. Among the rest, CO 2 efflux behaved either differentially from or similarly as pCO 2 depending on how k varied to Q in places of different stream hydraulics (Raymond et al. 2012). While simultaneous concentration and efflux decreases clearly indicate depletion of CO 2 sources as Q increased, presence of streams and rivers with decreasing pCO 2 but increasing CO 2 efflux suggests river hydraulics (and related surface-water turbulence status and k) was an equally important determining factor for CO 2 efflux-Q responses in rivers.
We examined the responsiveness of pCO 2 , k, and CO 2 efflux to Q along stream order (Fig. 3). Despite some regional variability (Fig. 3a), an increasing trend was identified for pCO 2 vs. Q with negative-to-neutral median slopes (−0.05 to 0) in small rivers of stream order 1-6 and significantly higher median slopes (0.02-0.24) in large rivers of stream order 7-10 ( Fig. 3a). Thus, stream order essentially captured the system dynamics where pCO 2 tended to decrease with Q in small streams and rivers but increase with Q in large rivers. Statistically, a positive median slope translated to more than half of positive responses at a specific stream order (vice versa). For example, positive pCO 2 -Q responses were 38-49% at stream order 1-6 and 56-100% at stream order 7-10, respectively (Supporting Information Table S2). We note large variations of observed slopes at each stream order. Mean slopes yielded a similar trend along stream order (Supporting Information Fig. S3a).
The responsiveness of k to Q was strongest in the smallest streams (stream order 1-3, with median slopes of~0.19), which decreased along stream order (Fig. 3b). The responsiveness of k was small at stream order 9-10 (with median slopes of~0.01) suggesting limited increase of k with Q in the largest river channels, although it should be noted that wind might play a role in large rivers with significant fetch. High responsiveness in small streams was likely due to steep stream slopes that favored quick dissipation of stream energy and strong increase of stream turbulence (Raymond et al. 2012). The observed decline in k-Q responses was in line with a decreasing k along stream order previously reported in US streams and rivers (Butman and Raymond 2011), suggesting stream slope not only exerted controls on k but also its variability in response to Q.
The responsiveness of CO 2 efflux to Q was positive (with median slopes of 0.12-0.33) across all stream orders. In contrast to the relatively unidirectional changes (either increase or decrease) of pCO 2 and k slopes along stream order (Fig. 3a,  b), the responsiveness of CO 2 efflux to Q was strongest at stream order 1-3 (with median slopes of 0.23-0.31), decreased along stream order 4-6 (with median slopes of 0.12-0.17), and increased again at stream order 7-10 (with median slopes of 0.13-0.33) (Fig. 3c). Strong increases of CO 2 efflux at stream order 1-3 were in sharp contrast with diluting CO 2 concentrations in these streams suggesting highly responsive k played a more decisive role in determining the CO 2 efflux-Q responses in these steeper terrains. Instead of transporting CO 2 downstream in their channels, these streams tended to emit a disproportionate amount of CO 2 to the atmosphere at higher Q. An examination across stream order indicates that these small streams were also associated with higher portions of opposite concentration and efflux responses than the other rivers (27-48% at stream order 1-4 vs. 0-26% at the other stream orders) (Supporting Information Table S2). Strong CO 2 efflux responsiveness at stream order 1-3 translated to high CO 2 48.7%

Liu and Raymond
Hydrologic controls on river CO 2 efflux input to these small streams considering high k and its responsiveness to Q would only lead to high CO 2 efflux if there was a large CO 2 source. In contrast, increased CO 2 efflux responsiveness in large rivers was related to stronger concentration increase considering limited increase in k with Q. Finally, weak increase in both CO 2 concentration and efflux at stream order 4-6 suggested reduced k and probably a smaller source supply to these rivers. We modeled how pCO 2 , k, and CO 2 efflux changed with Q (0.5-5 cm d −1 ) at different stream orders using the obtained log-linear relationships of each site (Fig. 4, median responsiveness shown). pCO 2 decreased with Q in streams and rivers of stream order 1-6 and increased in large rivers of stream order 7-10 ( Fig. 4a), in consistence with the slope observations (Fig. 3a). k increased with Q at all stream orders but most strongly in the smallest rivers (stream order 1-4) (Fig. 4b). CO 2 efflux increased with Q at all stream orders (Fig. 4c). Modeled CO 2 efflux was highest in small streams and decreased along stream order 1-6, suggesting decreasing terrestrial connectivity and CO 2 input (see discussion below). Modeled mean CO 2 efflux exhibited similar trends (see Supporting  Information Fig. S4).
A prior analysis (Hotchkiss et al. 2015) argued for a significant difference in CO 2 source supplies between small and large rivers in the US river networks. While terrestrial connectivity and lateral transport of landscape CO 2 comprised the largest CO 2 input to small streams and rivers, internal metabolism played a more important role in large rivers (Hotchkiss et al. 2015). Our analysis supports this presumption in that the responsiveness of CO 2 efflux to Q was strong at both ends of the river network while weak in medium-sized rivers (stream order 4-6) where neither of the two endpoint sources was prominent. While the continuous decreases in CO 2 efflux (Fig. 4c) and its Q responsiveness (Fig. 3c) along stream order 1-6 were in line with decreasing landscape connectivity and terrestrial CO 2 input in these rivers (Hotchkiss et al. 2015), internal river organic matter metabolism was still not a dominant process (Gomez-Gener et al. 2016;Winterdahl et al. 2016). Short-residence times and strong hydrologic responsiveness in these small rivers worked to shunt river organic matter quickly through the stream network (Raymond et al. 2016), leaving little opportunity for channel metabolism and processing until the higher orders were reached.
Our analysis further suggests that the "large river" effect where internal processes become more important (Hotchkiss et al. 2015) starts most likely at around stream order 6 or 7 where pCO 2 started to increase with and CO 2 efflux responded more actively to Q (Fig. 3). As internal processes played an increasingly important role in large rivers  (Hotchkiss et al. 2015), we propose strong responsiveness of pCO 2 and CO 2 efflux in these large rivers (Fig. 3a,c) was partly due to increased riverine metabolism at higher Q. Specifically, previous research (Raymond and Saiers 2010;Wilson et al. 2013) suggests increases of river organic matter concentration and especially its labile fractions at higher Q. Despite simultaneous changes in environmental conditions brought by shifted Q, increased substrate concentration would most likely facilitate aquatic organic matter metabolism and internal production of CO 2 (Lapierre et al. 2013;Rovelli et al. 2017). Most importantly, longer hydrologic residence times in these large rivers allowed for sufficient metabolic processing and longitudinal build-up of CO 2 in river channels (Raymond et al. 2016). Last, CO 2 uptake by primary production in the large rivers was limit due to high sediment loads at high Q (Dodds et al. 2013). Considering significant expansion of river corridor in large rivers (Thorp et al. 2006;McCluney et al. 2014), river-corridor interactions (which are an inherent attribute of Q dynamics in rivers) were expected to exert more controls on CO 2 -Q dynamics in large rivers than in small rivers. For example, more active hyporheic exchanges (Findlay 1995), reinundation of dried river channels, and re-connection of isolated riparian floodplains (Junk et al. 1989;Sieczko et al. 2016) are thought to supplement more respiration products to river channels under increased Q. However, to validate and systematically quantify this effect necessitate further research.
We do not rule out influence of the relative partitioning between CO 2 surface evasion and downstream transport, which was affected by k, on the hydrologic responses. Considering dominant CO 2 evasion over downstream transport as suggested by rivers of all sizes (or spatial scales) (Butman and Raymond 2011;Wallin et al. 2013), a small difference in the responsiveness of k to Q would have a major effect on pCO 2 or CO 2 efflux's responsiveness to Q. As an example, the responsiveness of k to Q in the largest rivers (stream order 9-10) was one twentieth of that in small streams (stream order 1-3) (Fig. 3). Thus, it is likely that the pCO 2 increase with Q in large rivers was (or partly) due to increased downstream transport under constrained emission. The same rule applies to rivers where k was highly responsive to Q (stream order 1-6) and rivers evaded a disproportionate amount of CO 2 under high Q, leading to negative to neutral concentration responses in these rivers. We propose a quantitative approach that takes into account both the concentration (downstream transport in channels) and evasion terms to quantify this effect.
We evaluated how pCO 2 and CO 2 efflux were distributed along an annual hydrograph by coupling the pCO 2 -Q and efflux-Q relationships with a daily Q dataset of 674 US Q Q Q k Fig. 4. Modeled changes of (a) pCO 2 , (b) k, and (c) CO 2 efflux with Q (0.5-5 cm d −1 ) at different stream orders. pCO 2 , k, and CO 2 efflux were modeled using corresponding log-linear relationships with Q at each site and median responses are shown here. stream/river sites (making up 82.8% of water quality dataset sites). We used daily Q observations of 10 yr (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) to reduce biases out of inter-annual Q variability. We found that flows greater than 2.3 cm d −1 (the median Q) were found in 9% of the days in a year, yet were responsible for 15% and 26% of the annual pCO 2 and CO 2 efflux, respectively (Fig. 5). Considering total CO 2 efflux from rivers was further determined by river surface area, greater contributions from high flows were expected when expansion of water surface area with Q was further taken into account. Future research is suggested to incorporate the hydrologic dynamics illustrated here for better modeling of CO 2 evasion from streams and rivers.

Conclusion
About half of conterminous US streams and rivers exhibited increasing pCO 2 -Q relationships and a significantly higher portion (80%) exhibited increasing CO 2 efflux-Q relationships. We found that small rivers (stream order 1-6) tended to evade rather than transport CO 2 downstream in response to stronger CO 2 inputs from terrestrial sources under increased Q. In large rivers, however both pCO 2 and CO 2 efflux increased with Q. We suggest that whereas decreasing CO 2 efflux responsiveness along stream order 1-6 was in line with a decreasing terrestrial input, positive pCO 2 and CO 2 efflux responses in large rivers suggested likely increased in-channel metabolism and strengthened rivercorridor exchanges at higher Q. Systematic changes in pCO 2 and CO 2 efflux's responsiveness to Q suggest surface evasion and downstream transport of CO 2 by river channels were closely related to stream order and hydrologic controls on both lateral terrestrial input and internal metabolism. We suggest to incorporate the hydrologic dynamics illustrated here for better constraining CO 2 emission from streams and rivers in future. Q over the annual hydrograph. The pCO 2 and CO 2 efflux fractions were calculated as a product between the relative frequencies of Q and modeled pCO 2 or CO 2 efflux values.