Landscape process domains drive patterns of CO2 evasion from river networks

Streams are important emitters of CO2 but extreme spatial variability in their physical properties can make upscaling very uncertain. Here, we determined critical drivers of stream CO2 evasion at scales from 30 to 400 m across a 52.5 km2 catchment in northern Sweden. We found that turbulent reaches never have elevated CO2 concentrations, while less turbulent locations can potentially support a broad range of CO2 concentrations, consistent with global observations. The predictability of stream pCO2 is greatly improved when we include a proxy for soil‐stream connectivity. Catchment topography shapes network patterns of evasion by creating hydrologically linked “domains” characterized by high water‐atmosphere exchange and/or strong soil‐stream connection. This template generates spatial variability in the drivers of CO2 evasion that can strongly bias regional and global estimates. To overcome this complexity, we provide the foundations of a mechanistic framework of CO2 evasion by considering how landscape process domains regulate transfer and supply.

CO 2 evasion is the largest pathway for carbon (C) loss from inland waters (Cole et al. 2007), with a total flux that matches net oceanic uptake. Particularly important are streams and rivers, which account for more than 85% of CO 2 evasion while only covering 17% of the inland water surface (Raymond et al. 2013). Yet, current global estimates of CO 2 evasion from running waters are remarkably unconstrained (Drake et al. 2017), with a threefold difference among studies, ranging from 0.7 (Lauerwald et al. 2015) to 1.8 Pg C yr −1 (Raymond et al. 2013), despite using the same data source (Hartmann et al. 2014). We suggest that this uncertainty is related to how we conceptually treat the two fundamental variables underlying CO 2 evasion: the partial pressure of water CO 2 (pCO 2 ) and the gas transfer velocity, k 600 , a measure that reflects water's turbulence.
There have been major advances in upscaling k 600 from fine scale measurements of stream morphology (Raymond et al. 2012), together with estimates of river area based on remote sensing (Allen and Pavelsky 2018). By comparison, we are still unable to predict pCO 2 and explain its variability over space and time. For instance, variability in pCO 2 is larger within small watersheds than across biomes (Wallin et al. 2018), and river size is a poor predictor of pCO 2 in the current global data set (Marx et al. 2017). To overcome this problem, current approaches for upscaling (Butman and Raymond 2011;Raymond et al. 2013;Borges et al. 2015;Wallin et al. 2018) typically involve estimating pCO 2 from aggregated water chemistry data over large regions and matching this with more finely resolved estimates of k 600 . However, summarizing pCO 2 and k 600 at different scales may result in erroneous estimates of CO 2 evasion because these underlying parameters may not be independent. Therefore, to upscale CO 2 evasion from streams, we need to develop robust tools to predict pCO 2 and k 600 at similar spatial scales.
Concentrations of CO 2 in streams reflect the dynamic balance of turbulent forces that promote evasion against multiple processes that resupply the system. On one hand, k 600 can vary greatly within meters depending whether it is estimated within a pool, riffle, or waterfall (Keller and Melhorn 1978;Wallin et al. 2011) and this will constrain the rate at which CO 2 degassing can occur; henceforth defined as the transfer limitation. On the other hand, CO 2 is continuously added to streams through groundwater inputs (Öquist et al. 2009) and internal production (Hotchkiss et al. 2015), which can also vary markedly at small temporal (Peter et al. 2014;Crawford et al. 2016) and spatial scales (Natchimuthu et al. 2017;Wallin et al. 2018;Lupon et al. 2019). This rate of CO 2 supply could thus also limit evasion to atmosphere; henceforth defined as supply limitation. The extent to which transfer vs. supply limitation regulate patterns of CO 2 evasion remains ill-defined, despite being fundamental to a mechanistic understanding of evasion across stream networks. Here, we ask: to what extent does transfer vs. supply limitation shape network-scale patterns of CO 2 evasion, and how does resolving these underlying drivers influence our ability to upscale evasion rates to larger areas?
To answer this question, we estimated CO 2 evasion and its drivers at reach scales for a stream network in northern Sweden (the 52.5 km 2 Miellajokka catchment, Fig. 1a). The catchment comprises a wide range of hydrological conditions and landscape properties and is thus suitable as a model system. We measured key hydrological parameters to estimate k 600 , directly determined the pCO 2 , and estimated CO 2 evasion across 168 sites in the stream network. Together with a quantification of wet areas in the catchment and other topographic information that might influence supply and/or transfer limitation, we assessed the drivers of the spatial variability in CO 2 evasion. We also utilize a global stream database (GLORICH; Hartmann et al. 2014) to test if our catchment-scale findings corroborate larger scale patterns observed for pCO 2 , k 600 , and evasion.

Site description
The Miellajokka catchment is in the arctic region in northwestern Sweden (68 21 0 14 00 N, 18 56 0 16 00 E). The catchment covers 52.5 km 2 , and elevation ranges between 384 and 1731 m above sea level (a.s.l.). The morphology of the catchment reflects several glaciations across millennia, which left behind U-shaped valleys and several deposits of coarse material and moraines. The broader landform of the catchment is characterized by steep peaks in the south, a large central plateau, and moderately sloped terrain that feeds an alluvial valley in the north (Fig. 1). The catchment is defined by four large streams which flow northward. The total length of the stream network is 44.6 km and includes first to fourth Strahler order streams (Supporting Information Table S1). Annual precipitation in Abisko (7 km from the study site, www. polar.se/abisko) for the period 1990-2013 was 337 mm, and the mean annual air temperature was 0.3 C. Land cover changes with elevation, from mostly bare rock above 1200 m a.s.l., to tundra vegetation between 700 and 1200 m, to sparse cover by mountain birch (Betula pubescens spp. Czerepanovii) and treeless tundra heath below the treeline at 700 m. Below 400 m, there are areas with more productive birch forest and denser canopy cover. Soil permafrost in this region is discontinuous at higher elevations and sporadic at low elevation areas. Bedrock is dominated by weakly metamorphic sedimentary rocks such as schist and dolomite in the central and western areas of the catchment, whereas more quarzitic and phyllitic hard schist and salic igneous rocks form bedrock along the eastern tributary.

Field sampling and analyses
We sampled 168 locations in this stream network between 11 th and 20 th of July 2016 (Fig. 1). The distance between sampling points increased with stream order as the travel time and velocity increases with channel size (Raymond et al. 2012); sampling locations were also constrained by having a suitable location to measure discharge (e.g., single channel, no pools). Thus, the distance between sites was between 30 and 100 m for first-order streams, between 100 and 200 m for second-order streams, and between 200 and 400 m for third-and fourth-order streams. Sampling was carried out between 08:30 h and 18:00 h every day. Weather conditions and flow levels across days were stable during this period, Rocher-Ros et al.
Landscape domains of CO 2 evasion from streams with the coefficient of variation (as percentage) of discharge at the catchment outlet of 12.2%. In addition, daytime variation in pCO 2 during sampling hours was monitored at seven sites and was, in general, constant during this time of day (Supporting Information Fig. S3). There was some variability between days but this was small compared to the spatial variability (Supporting Information Fig. S3). At each sampling site, we measured water depth and velocity every 0.2 m along transects using an electromagnetic flow meter (model 801 EC Meter; Valeport, Devon, U.K.). We also measured pCO 2 in situ using a handheld device (Vaisala DM70, Helsinki, Finland) adapted for wet environments as in Johnson et al. (2010) and water temperature with a handheld thermometer. Stream discharge was calculated as the product of the cross-section area and the velocity of the stream. CO 2 concentrations measured with the Vaisala instrument (accuracy of AE 10 ppm) need to be corrected for air pressure and temperature (Johnson et al. 2010). Specifically, the sensor output has to be corrected by 0.15% per hPa departure from 1013 hPa, and by −0.3% per C difference from 25 C, as stated in Johnson et al. (2010). The flux of CO 2 to the atmosphere was calculated as: where k CO2 is the gas transfer velocity for CO 2 , and ΔCO 2 is the difference in CO 2 concentration between the water and the air. We used an air CO 2 concentration of 380 ppm based on measured concentrations with the handheld device at the six monitoring sites. The air concentration was the same (380 ppm) at each site. The k CO 2 was calculated as: where Sc CO 2 is the Schmidt number for CO 2 at the water temperature calculated following the coefficients summarized in Raymond et al. (2012) and k 600 is the gas transfer velocity standardized for a Schmidt number of 600. Due to the impossibility of Rocher-Ros et al. Landscape domains of CO 2 evasion from streams directly measuring k 600 at a large spatial extent, we used measured hydraulic parameters to estimate the gas transfer velocity for each reach, based on the relationship following Raymond et al. (2012): where V is the water velocity (m s −1 ), S is the slope (unitless), Q is discharge (m 3 s −1 ), and D is depth (m). To obtain k 600 for the global database GLORICH, we derived the slope for each site and calculated k 600 using velocity and slope as predictors. The pCO 2 used from the GLORICH database (Hartmann et al. 2014) was processed as follows: first we calculated the median value if several per site were available, then we excluded sites with unclear geographical coordinates and/or with pCO 2 > 40,000 ppm, which gave us a total of 4790 sites. A full description of the data used as well as the modeling of k 600 can be found in the Supporting Information Text S1. We used the digital elevation model at a resolution of 2 m to derive a depth to water index for the Miellajokka Catchment (DTW; details of the computations can be found in Supporting Information Text S2). Wet areas were defined as locations where the distance to the water table was less than 1 m, suggesting stronger connectivity to the stream network (Lidberg et al. 2017). For each site, wet area was computed as the percent cover of the incremental catchment area from the previous upstream site with a DTW < 1 m. A source of error when estimating the effect of wet areas for a given reach is that the travel distance of a gas may be longer than the distance to the upstream sampling site, and therefore integrate the properties from that upstream site. We calculated the mean travel distance of a gas entering the stream as 0.7 × V/K 600 (V in m d −1 , and K 600 is k 600 /D in d −1 ) following Demars et al. (2015), and found that the median integration distances for first-, second-, third-, and fourth-order streams were 44 m, 78 m, 86 m, and 678 m, respectively. These values are similar or smaller than the distance between the sampling sites.

Statistical analysis
We completed all statistical analyses in R (version 3.5.1), using the lm() function for linear and multiple regressions, and ggplot2 (version 3.0.0) for visualization. To assess the relative importance of k 600 and pCO 2 on CO 2 evasion, we used a multiple linear regression with the variables standardized and log-transformed. Accordingly, Eq. 1 defined in log-space becomes: log F CO 2 ð Þ= log k CO 2 ð Þ+ log ΔCO 2 ð Þ , and in a multiple linear regression (i.e., log , the coefficients x 1 and x 2 are a measure of the relative importance of k CO 2 and ΔCO 2 to CO 2 evasion. To do this, we first standardized variables using the scale() function in R. Then, to calculate the partial R 2 for each independent variable, we used the function modelEffectSizes() of the package lmSupport (version 2.9.13), which provides the semi-partial correlation coefficients (r) and the partial coefficients of determination (R 2 ) for the two independent variables. Finally, we used a bootstrapping exercise to explore the sampling intensity required to capture median CO 2 evasion from Miellajokka streams. Here, we randomly sampled the population of measurements assuming a variable sample size, from 1 to 168 (with sample replacement). This sampling was repeated 200 times for each sample size, and for each simulation, we calculated the median CO 2 evasion rate. The data sets used in the study are deposited in the Swedish National Data Service (https://doi.org/10.5878/pxa2-vy55) (Rocher-Ros 2019), and code to reproduce the results and figures is deposited in GitHub (https://github.com/rocher-ros/co2_domains_publication).

Controls on CO 2 evasion from streams
In a worldwide data set of stream water chemistry (Hartmann et al. 2014), we find that CO 2 concentrations are always low (< 1000 ppm) at k 600 > 50 m d −1 , while there is large variation in stream pCO 2 below a k 600 of 50 m d −1 (Fig. 2a). Thus, turbulent streams never have higher CO 2 concentrations, while less turbulent streams can potentially support a broad range of concentrations, from being in equilibrium with the atmosphere to concentrations above 20,000 ppm. Our empirical estimates from the Miellajokka catchment are remarkably consistent with this global observation. Across the 45 km long stream network (Fig. 1a,b), we consistently observed low CO 2 concentrations (< 500 ppm) at more turbulent stream reaches (k 600 > 100 m d −1 ) while high pCO 2 values (> 1000 ppm) only were observed at a k 600 below 100 m d −1 (Fig. 3a). The absence of high pCO 2 values at the more turbulent stream reaches suggests that the degassing rate is considerably larger than the supply of CO 2 that supports the evasion to the atmosphere. Hence, evasion from such stream reaches is likely to be supply rather than transfer limited as has been shown in other mountainous areas (Crawford et al. 2015).
While supply limitation prevents the build-up of higher stream CO 2 concentration in turbulent reaches, we also found that where k 600 < 100 m d −1 stream pCO 2 exhibited large variation, with both high and low values (Fig. 3a). Thus, at similar transfer rates (e.g., similar k 600 ), the supply of CO 2 to streams varied markedly across this network. To explore the controls over this variability, we used the percentage of wet area draining into each site as a proxy for the strength of stream-groundwater connectivity (Kuglerová et al. 2014;Lidberg et al. 2017), assuming this also reflects the supply of C from riparian zones to streams (Hotchkiss et al. 2015;Tank et al. 2018;Lupon et al. 2019). We found that below a threshold of k 600 = 42 m d −1 (Supporting Information Fig. S5), stream pCO 2 was significantly positively related to the % wet area in the adjacent riparian zone and this explained 63% of the variation in CO 2 concentrations (Fig. 3b). This relationship suggests a strong component of supply limitation on CO 2 evasion-even in streams with low turbulence, where CO 2 evasion should be primarily transfer limited. In contrast, only 23% of the variation in pCO 2 could be explained by the % wet area at k 600 values above 42 m d −1 (Fig. 3b). This second relationship most likely reflects the overarching influence of degassing, which at high exchange rates is able to mask the effects of variable supply to streams.
The balance between the supply and transfer shapes the overall relationship between k 600 and pCO 2 (Fig. 3a). The nature of this relationship limits the predictability of pCO 2 when solely  based on k 600 even if the two parameters are related (Fig. 3a). Indeed, including % wet area in a multiple linear regression greatly increased our ability to predict pCO 2 across all sites from 32% to 62% (pCO 2 (ppm) = 470 − 1.7 × k 600 + 17.9 × % wet area, p < 0.001). While we lack data to test this globally at such a fine scale, both catchment ) and continental scale (Borges et al. 2015) studies have similarly shown that wetland cover can explain patterns in riverine CO 2 . Regardless, it is clear that the relationship we see between pCO 2 and k 600 within a 52 km 2 catchment also holds at a global scale (Fig. 2a), suggesting a universal threshold of a k 600 beyond which turbulent flows prevent the build-up of CO 2 . Understanding whether pCO 2 or k 600 drives CO 2 evasion is necessary to upscale CO 2 emissions globally. Surprisingly, in the Miellajokka catchment, hotspots for CO 2 evasion are primarily found in stream reaches with lower k 600 and not in the more turbulent reaches of the catchment (Fig. 3c). Importantly, CO 2 evasion rates are strongly related to pCO 2 in this catchment (Fig. 3d) and in the global data set (Fig. 2c). The importance of pCO 2 for evasion is also clear in the multiple linear regression. For the Miellajokka, the multiple linear regression was log F CO 2 ð Þ= 0:62 × log k CO 2 ð Þ+ 1:27 × log ΔCO 2 ð Þ (R 2 = 0.98, p < 0.001), with a partial R 2 of 0.89 for log(ΔCO 2 ) and 0.04 for log k CO2 ð Þ. Thus, the higher coefficient (1.27 vs. 0.62) and the higher R 2 (0.89 vs. 0.04) of log(ΔCO 2 ) compared to log k CO 2 ð Þ imply that pCO 2 is the critical driver of evasion and not k 600 . Results from global data set were similar; here, the multiple linear regression was log F CO 2 ð Þ= 0:42 × log k CO 2 ð Þ+ 1:01 × log ΔCO 2 ð Þ(R 2 = 0.97, p < 0.001). The partial R 2 for log(ΔCO 2 ) was 0.89 and for log k CO 2 ð Þwas 0.02. Together, these results highlight that constraining pCO 2 at fine spatial scales is essential to determine and upscale CO 2 evasion rates.
A turbulent stream reach is thus not a hotspot for CO 2 evasion by itself; our data suggest that such locations are instead found at reaches with moderate to low turbulence, but where there is sufficient supply of CO 2 to maintain higher stream concentrations and thus atmospheric losses. Importantly, while the supply of CO 2 to the water column emerges as the main driver of CO 2 evasion in these streams, the processes behind its origin are still uncertain (Hotchkiss et al. 2015). Some studies show that direct inputs of CO 2 produced in soils dominate this flux (Öquist et al. 2009), others highlight that near stream, riparian soils are the major source (Leith et al. 2015), while others still suggest that instream aquatic respiration can account for a large proportion of the CO 2 produced (Rasilo et al. 2016;Demars 2018). These are not mutually exclusive processes, and we are not able to resolve among them, as our measure of wetness index only captures variation in hydrological connectivity between land and water and changes in the lateral extent of riparian soils. These wet areas may foster higher rates of soil metabolism (Riveros-Iregui and McGlynn 2009) and directly supply streams with CO 2 . At the same time, these zones may yield higher rates of lateral organic C supply, and thus support higher rates of hyporheic (Pusch 1996) and/or benthic (Battin et al. 2008) respiration.

Process domains and CO 2 evasion in river networks
Overall, our results show that transfer limitation only regulates where CO 2 evades, while supply ultimately controls the magnitude of evasion. In the Miellajokka catchment, reaches with steeper channel slopes have greater turbulence and thus potentially elevated evasion rates, but only if there is sufficient CO 2 from internal production (e.g., aquatic respiration) or from lateral or upstream sources. By contrast, in many lower gradient reaches, riparian and groundwater features promote greater lateral C flux to streams and potentially elevated rates of aquatic metabolism and CO 2 production. Yet, such zones often lack the kinetic energy to support very high rates of CO 2 evasion. Thus, the true hotspots for evasion in the network likely emerge at the transitions between reaches characterized by high C supply and those by high turbulence-where high values of pCO 2 and k 600 may coexist (Fig. 4). This perspective places new emphasis on longitudinal connectivity among transfer-and supply-limited reaches as a driver of C gas efflux. It has been suggested that the connectivity between lotic and lentic environments in the aquatic continuum creates a similar pattern (Hotchkiss et al. 2018), but here we suggest that these spatial interactions may operate at even finer scales within river networks.
While these spatial interactions likely generate unique patterns of CO 2 evasion for every catchment, the broader role of topography as the overarching driver of supply vs. transfer limitation in streams and rivers is general (Fig. 4). It is well recognized that shifts in topography along river systems generate discontinuities in geomorphic structure, or process domains that are characterized by unique hydrologic, edaphic, and ecological dynamics (Montgomery 1999). Moreover, recent studies have highlighted the fundamental roles that geomorphic structure plays as a regulator of C inputs (Duvert et al. 2018;Lupon et al. 2019), storage (Wohl et al. 2012), biological processing (Jankowski et al. 2014), and transport in river systems (Battin et al. 2008). Here, we provide evidence that simple topographic tools can help distinguish domains in the landscape where lateral C supply to streams is the dominant process driving gas evasion from others where turbulent flow and degassing are paramount. The transitions in these dominant processes underpin the fundamental dependence between pCO 2 and k 600 that appears consistent at a global scale (Fig. 2a). Finally, we suggest that the relative position of these domains determines the location of hotspots for evasion in landscapes, which can also cause an extreme bias when upscaling gas evasion at regional or global scales. Our framework adds to the general understanding of C cycling in river networks, but also provides a set of guiding principles that will aid efforts to scale up CO 2 emissions from these systems.

Implications for upscaling
The effect of landscape domains can be very important when upscaling CO 2 evasion. Within the Miellajokka catchment CO 2 evasion ranged from 0 to 60 g C m 2 d −1 , the latter among the highest evasion rates measured in northern streams (Natchimuthu et al. 2017) (Fig. 5b). Importantly, the distribution of these rates was highly skewed toward lower values (Fig. 5b), resulting in a large difference between the median and the mean estimate, 3.34 g C m −2 d −1 and 7.32 g C m −2 d −1 , respectively. To what extent does this spatial heterogeneity in evasion rates also affect our ability to predict broader, catchment-scale CO 2 evasion estimates? The bootstrapping of the distribution of CO 2 evasion (Fig. 5a) shows that by collecting 10 samples from streams in this catchment would result in a large uncertainty, with a 5 th -95 th percentile of 1.43-8.02 g C m −2 d −1 . Collecting 50 samples narrows this uncertainty to a 5 th -95 th percentile of 2.23-4.86 g C m −2 d −1 , closer to the observed median of 3.34 g C m −2 d −1 . As the catchment covers 52.5 km 2 , these results then indicate that constraining catchment-scale estimates of CO 2 evasion would require roughly one sample for each square kilometer of catchment.
While the spatial variation in pCO 2 and k 600 is incorporated into the above estimates of CO 2 evasion, most current upscaling studies rely on estimates of pCO 2 , k 600 , and stream area that are determined at different scales. If we apply a similar approach to the Miellajokka catchment and use the median k 600 determined from stream morphology across the catchment and pCO 2 measured at the main stem outlet (which is also a national monitoring station), the estimated CO 2 evasion is 54.7 g C m −2 d −1 (red triangle in Fig. 5a), more than 15-fold higher than the median flux. These differences in the parameters to estimate CO 2 evasion will thus result in a grossly overestimated catchment CO 2 evasion. In our case, we overestimate evasion, due to the high median catchment k 600 (54.5 m d −1 ) and high pCO 2 in the monitoring station as it is located in a flat and moist area (1180 ppm, site M1 in Supporting Information Fig. S1), but the direction of this bias could go either way. This exercise highlights that the scale at which parameters are estimated can strongly affect the evasion estimate. For example, Raymond et al. (2012) calculated riverine CO 2 evasion globally using pCO 2 from the GLORICH database but used topographic data at finer spatial scales to determine the k 600 , obtaining a CO 2 evasion rate of 3358 g C m 2 yr −1 (total riverine CO 2 evasion divided by the reported riverine area). In this study, we calculated the k 600 and the CO 2 evasion at the sites of the GLORICH database where pCO 2 was measured, obtaining a median evasion rate of 1670 g C m 2 yr −1 , about half of the estimate by Raymond et al. (2012). Our measure is much closer to the estimate of riverine CO 2 evasion of 1574 g C m 2 yr −1 from Lauerwald et al. (2015), which upscaled pCO 2 and k 600 at the same spatial scales although Fig. 4. Landscape process domains shaping CO 2 concentration and evasion rates in streams. pCO 2 can respond strongly to inputs of C from land or to outputs via evasion. The spatial arrangement of C inputs and turbulent areas determines the controls and the magnitude of CO 2 evasion. Thus, CO 2 evasion will be either transfer limited (A), supply limited (C) or both (D). Hotspots of CO 2 evasion can then occur when supply is high and there is enough turbulence to degas the water (B).
at spatial resolution too low to capture the variations we described (about 3000 km 2 ). This indicates that the global estimates of riverine CO 2 evasion may be overestimated if they fail to capture spatial patterns of pCO 2 and k 600 at the same fine spatial scales.
In summary, we find three critical shortcomings that hinder our ability to generate accurate predictions of CO 2 evasion from streams: first, the high spatial variability per se in CO 2 evasion and its parameters across stream networks. Second, the use of different scales to summarize the main variables that predict CO 2 evasion, and third, the lack of tools to predict pCO 2 across stream networks. Thus, our results suggest that the mismatch between the current global estimates of riverine CO 2 evasion (0.7-1.8 Pg C yr −1 ; Raymond et al. 2013;Lauerwald et al. 2015), and the recent reviewed value of 3.9 Pg C yr −1 (Drake et al. 2017) is partly due to the shortcomings we identified. We are aware of the unfeasibility of an extensive sampling program at the spatial scale that we identified (1 km 2 ). Therefore, we suggest that, to advance our knowledge of the role of inland waters in the C cycle, we should move beyond accounting techniques of upscaling, and rather develop mechanistic models anchored by topographic features that govern both inputs and exchange of CO 2 (Fig. 4). The red triangle shows the CO 2 evasion if calculated with the measured pCO 2 at the monitoring station close to the outlet, and median k 600 of the catchment. The red dashed line shows the median CO 2 evasion, also in the histogram of the observed CO 2 evasion rates (b).