Aquatic‐Terrestrial Linkages Control Metabolism and Carbon Dynamics in a Mid‐Sized, Urban Stream Influenced by Snowmelt

Freshwater streams can exchange nutrients and carbon with the surrounding terrestrial environment through various mechanisms including physical erosion, flooding, leaf drop, and snowmelt. These aquatic‐terrestrial interactions are crucial in carbon mobilization, transformation, ecosystem productivity, and have important implications for the role of freshwater ecosystems in the global carbon budget. We utilized high‐frequency oxygen, temperature, and carbon dioxide (CO2) data to infer watershed connectivity in Boulder Creek, a mid‐sized (1160 km2) watershed located in Colorado, USA. Daily modeled gross primary production (GPP), ecosystem respiration (ER), net ecosystem production (NEP), and reaeration coefficients (K600) were paired with high‐frequency, in‐situ dissolved CO2 data to characterize changes in metabolic regime and carbon flux on a stream influenced by seasonal snowmelt. GPP and ER were correlated (ρ = −0.72, p ≪ 0.001) during the non‐snowmelt period and NEP was frequently negative. Mean FCO2 during the non‐snowmelt period was approximately 302 (±171) mmol C m−2 d−1 and was primarily supported by watershed CO2 inputs. During snowmelt, GPP and ER were not significantly correlated (ρ = −0.22, p = 0.05), and mean NEP was significantly more negative than during non‐snowmelt. Watershed connectivity was higher during snowmelt, as evidenced by significantly higher FCO2 (843 ± 338 mmol C m−2 d−1) and greater allochthonous CO2 inputs than during non‐snowmelt periods, emphasizing the effects of seasonal differences in aquatic‐terrestrial linkages in this stream. We suggest that our understanding of watershed carbon budgets is subject to temporal dynamics which control the degree of connectivity between terrestrial and aquatic ecosystems.

The influence of terrestrial inputs on stream F CO2 is related to stream size and discharge, which creates spatial variability and influences the magnitude of F CO2 from streams (Hotchkiss et al., 2015). Smaller, low order, streams have relatively high loading of terrestrial carbon, both OC and CO 2 , which significantly influences the carbon and metabolic dynamics of the stream, resulting in a high rate of F CO2 . In contrast, the relative influence of terrestrial carbon on larger rivers, or high order streams, is diminished, and CO 2 dynamics are governed to a larger degree by processes internal to the stream, such as NEP (Hotchkiss et al., 2015). However, the spatial structure of stream carbon biogeochemistry and its relation to metabolism is likely to be temporally variable due to differences in the degree of watershed connectivity. Therefore, expanded understanding of temporal and spatial controls on stream F CO2 and its relation to internal and external controls on the watershed is needed, as it can ultimately impact water quality and the global carbon budget (Butman & Raymond, 2011;Butman et al., 2016).
New tools, especially relatively inexpensive and reliable environmental sensors and modeling approaches, allow for long-term and expanded insights into metabolic rates, carbon fluxes, temporal stream dynamics and drivers (Battin et al., 2009;Bernhardt et al., 2018;Smith & Kaushal, 2015). The expanding prevalence of high-frequency CO 2 observations allows more detailed analysis of F CO2 . Pairing these observations with dissolved oxygen (DO) also allows the estimation of the gas exchange velocity from metabolism model estimates. These estimates can be made on a daily timescale, which contrasts with approaches that rely on periodic direct measurements of F CO2 from stream surfaces (Striegl et al., 2012) or broad syntheses that couple calculated pCO 2 with statistical representations of the gas exchange velocity (Battin et al., 2009;Butman & Raymond, 2011). Current F CO2 estimates from fresh waterbodies are based on relatively sparse data and typically do not account for diel changes in dissolved CO 2 concentrations (Battin et al., 2008;Butman & Raymond, 2011;Butman et al., 2016;Stets et al., 2017). Therefore, a greater understanding of the effects of seasonal dynamic hydrology on terrestrial productivity and aquatic ecosystem processes can be gained by incorporating updated CO 2 observations and modeling techniques to include high-frequency F CO2 estimates.
We used high-frequency DO and CO 2 observations to gain expanded insights into the seasonal patterns of stream metabolism and F CO2 and the implications for watershed connectivity in a snowmelt-dominated stream. We sought to answer the questions: How does the relationship between GPP and ER, including the rate of NEP, vary between snowmelt and nonsnowmelt seasons? Does the magnitude of F CO2 vary seasonally? How much of the F CO2 can be attributed to stream net heterotrophy (negative NEP) and does this vary seasonally? High-frequency sensor data, including DO, pCO 2 , light, and temperature, coupled with stream discharge data, were collected over a 2-year period in Boulder Creek, Colorado, USA. Daily modeled metabolism rates and high-frequency F CO2 estimates were crucial to identify seasonal changes in stream productivity and ultimately carbon transformation in the stream.

Site Overview
The Boulder Creek watershed is approximately 1,160 km 2 and is located on the eastern slope of the Rocky Mountains in Colorado (Figure 1). The headwaters originate at the continental divide and flow through the alpine, subalpine, montane, foothills, and plains ecoregions. The upstream, mountainous part of the watershed is sparsely populated and primarily forested, whereas the lower plains region consists of urban and agricultural land use (Murphy, 2006). The study site was co-located with U.S. Geological Survey (USGS) stream-gaging station 06730200 (USGS, 2020), on Boulder Creek downstream of the city of Boulder (but upstream of the city's wastewater treatment plant), at an elevation of 1,562 m above sea level with a 795 km 2 drainage area. Only 8% of the Boulder Creek Watershed upstream of the study site is developed (urban) land use, indicated with gray shading in Figure 1 (calculated using StreamStats). The area near the monitoring location is used recreationally for fishing, swimming, biking, and walking.

Discharge and Stream Depth
Seasonal and annual variations in stream discharge (Q) of Boulder Creek are dominated by snowmelt events and water management practices (Murphy et al., 2003). Snowmelt events are defined by seasonal increases in Q caused by rapidly melting snow in the warmer spring and summer months, and therefore were determined from the rate and timing of streamflow of Boulder Creek ( Figure S1). The beginning of snowmelt is evident from a rapid, continuous rise in Q in the spring. The end of the snowmelt period was estimated from timing of a rapid decrease to a relatively constant Q at an upstream gaging station with fewer diversions (Colorado Division of Water Resources State Gage BOCOROCO near Orodell, CO). During the snowmelt season, melting alpine snowpack exceeds the holding capacity of reservoirs in the basin and exceeds the rate of diversions in the watershed (Murphy, 2006), which allows a stronger connection between the upper and lower parts of the basin. The snowmelt period (approximately April to July) is also the wettest time of the year in the lower part of the basin. For example, on average, half of the 50 cm of precipitation that Boulder receives annually occurs during the period between April and July (https://www.esrl.noaa. gov/psd/boulder/Boulder.mm.precip.html), as snow, rain, or mixed precipitation. Water management practices that alter the flow regime of Boulder Creek include reservoirs, diversions into and out of the watershed, and wastewater effluent inputs. A trans-basin supply canal located about 800 m upstream of the study site (referred to here as the "supply canal"; Figure 1) delivers water diverted from the western slope of the Rocky Mountains to Boulder Creek from April to October (Murphy, 2006), and contributed between 0% and 95% of the base flow in Boulder Creek during the study period ( Figure S2). Discharge data for the supply canal (Northern Water, https://www.northernwater.org) minus a diversion from the canal was used to determine actual streamflow contributions to Boulder Creek. This input could have affected the DO and dissolved CO 2 dynamics at the study location; therefore, we limited the metabolism analysis to days when this diversion contributed less than 25% of the total flow of Boulder Creek ( Figure S2). Overall, 110 days of DO data and 88 days of pCO 2 data were removed from the analysis due to >25% of Boulder Creek streamflow originating from the canal.
Stream discharge measurements for Boulder Creek at 75th Street (USGS stream-gaging station 06730200) were obtained from the USGS National Water Information System (USGS, 2020). Mean reach depth was measured in November of 2016 and was related to the gage height, which is reported at 15-min intervals. Mean reach depth was estimated from cross-sections performed every 10 m over the calculated DO reach length, determined as 3v/K 600 , where K 600 is the reaeration coefficient and v is average velocity in m d −1 (Chapra & DiToro, 1991). After several months of data collection in which initial estimates of GPP, ER, and K 600 were obtained (not shown), the model was re-run for the entire data collection period after correcting to the new mean depth estimate.

Dissolved Gases, Temperature, and Light Measurements
Temperature, DO, dissolved CO 2 , and light intensity measurements were collected from Boulder Creek upstream of 75th Street from June 2016 to May 2018. Temperature and DO measurements were collected at 10-min intervals using a submersible miniDOT dissolved oxygen logger manufactured by Precision Measurement Engineering (Vista, CA, USA). The factory calibration for the miniDOT probes was internally validated and subsequently utilized for the first year of data collection. Precision Measurement Engineering ensured DO and temperature calibration to be within ±5% and 0.1°C, respectively. The second year of use, the miniDOT probes were calibrated using a membrane inlet mass spectrometer (MIMS; Hall et al., 2016). Calibration samples were collected by syphoning Boulder Creek water into a 12-mL glass exetainer in triplicate and preserving with 100 μL of 50% w/v zinc chloride, then carefully capping to ensure no air bubbles were present in the vial. Samples were stored at room temperature until analysis. Correction factors were calculated by converting the mass-32 peak obtained from the MIMS software into a DO concentration and comparing the value to the DO measurement in the stream when the sample was collected (adjusting for temperature and pressure; Kana et al., 1994). Calibration factors were assumed to be a linear offset and were between 0.06 and 0.73 mg L −1 in year two of the study. The temperature sensor was factory calibrated and was assumed to maintain an error range of ±0.1°C throughout the experiment.
Dissolved CO 2 measurements were collected using a Vaisala GMT220 with a CO 2 concentration range of 0-10,000 ppm CO 2 (Vantaa, Finland). The probe was wrapped in high-porosity polytetrafluorethylene (PTFE) sleeve (a gas-permeable, water-impermeable material) sealed with plasticized rubber to prevent water entry (Johnson et al., 2010). The probe was calibrated with 0, 390, 1051, 2024, 4020, and 6080 ppm certified CO 2 standards purchased from Airgas. The calibration curve related a known CO 2 concentration to a voltage output. The slope and offset generated by the standard curve were incorporated into the SEVolt command of the Campbell Scientific CR1000 code (Crawford et al., 2017). Dissolved CO 2 measurements were recorded using the Campbell Scientific CR1000 Measurement and Controller System with 1-h resolution. Raw CO 2 concentrations were adjusted for temperature and pressure and converted to mole fraction CO 2 . The partial pressure of CO 2 in the aqueous phase (pCO 2 ) was calculated by multiplying the mole fraction by the barometric pressure in atmospheres. Barometric pressure was estimated to be 850 hPa for the entire study based on the elevation at the study site. We have found in previous intercomparisons of several CO 2 sensing/measurement methods that the Vaisala GMT220 probe performed well when compared to serum bottle, syringe equilibration, or carbonate system estimates (calculated from alkalinity, temperature, and pH) under a range of conditions. Stets et al. (2017) reported strong coherence between Vaisala CO 2 measurements and those of the syringe equilibration method (In a comparison done for another project, we found: Vaisala CO 2 (ppm) = 576 + 0.98(Serum CO 2 ), rsq = 0.97 (Stets et al., 2017).
Light intensity measurements were collected using an Onset HOBO Pendant Temperature/Light Data Logger (Cape Cod, MA, USA). The pendants were factory calibrated at the start of the experiment. According to manufacturer specifications, light measurements are relative readings and account for light waves between 500 and 1100 nm with ≥50% response. Modeled light was used instead of light measurements if: (1) stream depth was ≥1.59 m; (2) if an observation between 10:00 and 15:00 local time registered as zero; or, (3) if there were no observed light data due to equipment malfunction. Approximately 50% of the light data utilized in the metabolism model were modeled light. Light data were fit to a localized polynomial regression to minimize outliers.
All probes were attached to a steel fireplace grate and securely positioned in the thalweg of the streambed with rebar. All sensors were fully submersed for the duration of the study and were oriented downstream to prevent damage from large floating objects. The light sensor was oriented horizontal to the streambed. Sensors were maintained on a weekly basis during warmer months (approximately April-September), and biweekly basis during colder months (approximately October-March). Gaps in the final dataset are due to either sensor malfunction or exclusion of the data due to the supply canal providing >25% of streamflow to Boulder Creek at the study site.

The Metabolism Model
Stream metabolism was estimated using the R package streamMetabolizer (Appling et al., 2017). The model uses a hierarchical Bayesian approach to solve the equation: where O 2t , P t , R t , and F O2t represent DO concentration and volumetric expressions of gross photosynthesis, ER, and the change in O 2 across the air-water interface at time t, respectively. R t is negative by extension of ER being designated as negative. The model fits the volumetric parameters GPP, ER, and the reaeration coefficient K O2 to reproduce the oxygen time series as closely as possible using an inverse modeling approach. R t was related to ER as where z represents the depth of the water column and R t and ER are volumetric units.
P t was related to GPP through z and photosynthetic photon flux density (PPFD) as where PPFD t is the photosynthetic photon flux density at time t, measured as light intensity, and ∑PPFD is the total light flux for a 24-h period. Because ∑PPFD is required to run the model, light data were handled such that days with complete observational records were used preferentially. On some days incomplete or unusable light data were collected due to instrument failure or extreme shading from sediment, detritus, or debris. In these cases, PPFD was calculated using geographic location and day of the year (Yard et al., 2005) and is a built-in function of streamMetabolizer (Appling et al., 2017).
We assumed a linear response of stream GPP to light in the metabolism model. Modeling experimentation showed that the linear light model was conceptually appropriate for a natural stream setting due to the bulk of metabolism occurring at the benthic surface, where photosynthetic saturation is less probable due to light stratification of the water column. We tested this assumption by running the model in both linear and saturated light configurations and superimposing the modeled and observed DO data. The linear light model yielded more precise DO estimates, especially on highly productive days where maximum DO concentrations exceeded 12 mg L −1 . In this system, the saturated light configuration underestimated the daily maximum dissolved oxygen concentration as much as 25% and in some extreme cases overestimated by 150%.
In addition to GPP and ER estimates, streamMetabolizer estimates the flux of O 2 across the air-water interface and is given by where K O2 is the oxygen reaeration coefficient, O 2A is the O 2 concentration at equilibrium with the atmosphere, and O 2t is the observed O 2 concentration at time t. K O2 was related to K 600 using the temperature-dependent Schmidt number (Jähne et al., 1987). O 2A was calculated for each time interval (t) based on temperature and barometric pressure.
Parameter values and uncertainties were obtained using Markov Chain Monte Carlo (MCMC) including both process and observation error, which helps to guard against overfitting (Appling et al., 2018). Partial pooling was achieved by binning daily discharge values to constrain K O2 . Prior estimates for GPP were set to 93 (±63) mmol C m −2 d −1 based on preparatory experiments estimating metabolic rates in Boulder Creek. Two thousand burn-in steps were run prior to model parameter value output collection. All model options were specified in streamMetabolizer using the specs function. Model estimates were converted from g O 2 m −2 d −1 to mmol C m −2 d −1 by assuming a one-to-one molar ratio between O 2 and C (Stets et al., 2017). Tests of seasonal difference were subject to autocorrelation bias because of the unequal sample density between seasons, so t-test values were corrected using Newey-West adjusted standard errors using the sandwich package in R (Newey & West, 1986).

Calculations and Estimates of F CO2
Values for pCO 2 were converted to concentration values using Henry's law (Butler, 1982). F CO2 at time t (F CO2t ) was also estimated using the calculated CO 2 concentrations and the K 600 values produced through streamMetabolizer using the equation where kCO 2 is the gas exchange velocity of CO 2 , CO 2A is the concentration of CO 2 in equilibrium with the atmosphere at time t, and CO 2t , is the observed CO 2 concentration at time t. kCO 2 was calculated from K 600 using stream depth (z) and the temperature-dependent Schmidt number (Jähne et al., 1987;Striegl et al., 2012). CO 2A was calculated from Henry's law constant by assuming that atmospheric CO 2 was 390 ppm and correcting for barometric pressure. Hourly F CO2t estimates were integrated to produce daily estimates of F CO2 . By convention, positive values of F CO2t indicate a mass loss of CO 2 from the stream to the atmosphere.
Pairing daily, modeled K 600 with pCO 2 sensor readings to estimate F CO2 minimizes artifacts associated with the carbonate buffering system, such as diel lags in CO 2 concentration and flux (Stets et al., 2017), and expanded the temporal scale such that F CO2 could be reasonably compared to daily NEP to determine the contribution of NEP to F CO2 in Boulder Creek.

Hydrologic Regime
The study period encompassed one full snowmelt period and two partial snowmelt periods (Figure 2). Mean Q during both of the snowmelt periods was approximately eight times higher than during the non-snowmelt periods, with over half of annual discharge (54%) occurring in May and June (Figure 2). Peak discharge in 2016, 2017, and 2018 was 23.8, 23.5, and 24.7 m 3 s −1 , respectively, compared to 34.3 and 41.3 m 3 s −1 in 2014 and 2015 respectively (Figure 2). A statistical summary of Q for the full study can be found in Table 1.

DO and pCO 2 Dynamics
The DO record covers 518 days out of the 723 total days of the study. Of these 518 days, 110 days (21%) were removed from the analysis due to >25% of streamflow originating from the supply canal and 53 days (10%) contained insufficient data for the metabolism model and were dropped from analysis. Overall, 355 days (49% of the study period) are represented in the final metabolism analysis (Figure 3a). The relationship between observed DO and pCO 2 can be found in Figure S3. The record for pCO 2 covers 209 days out of the 723 total days of the study (29%). However, we removed 88 days from the analysis due to (1) the supply canal contributing >25% of baseflow to Boulder Creek at the study site and (2) non-overlapping modeled metabolism rates and pCO 2 data. The final analysis contained 121 days of pCO 2 data ( Figure 3c).  Figure S1. Note. Values account for observations across the entire study. Diel patterns in pCO 2 persisted in both snowmelt and non-snowmelt conditions, with less prominent fluctuations observed in the spring and larger fluctuations in the fall (Figure 3d). Mean pCO 2 across the entire study was approximately 1,500 (±531) ppm ( Figure 3c). Monthly mean pCO 2 was highest in October (1,690 ± 680 ppm) and lowest in June (1,410 ± 465 ppm). Mean pCO 2 was slightly lower during snowmelt (1,480 ± 431 ppm) than in the non-snowmelt (1,580 ± 613 ppm) period. Carbon dioxide flux (F CO2 ) was also estimated as the product of observed CO 2 concentrations and K 600 obtained from the metabolism model. Mean F CO2 across the entire study was 506 (±363) mmol C m −2 d −1 (Figure 7b). Mean F CO2 during the snowmelt period (843 ± 338 mmol C m −2 d −1 ) was almost three times higher than mean F CO2 during the non-snowmelt period (302 ± 171 mmol C m −2 d −1 ).

Metabolism Model Estimates: GPP, ER, NEP, and K 600
Across both years of the study, January had the lowest mean rate in ER (157 mmol C m −2 d −1 ), while May had the highest mean rate in ER (−458 mmol C m −2 d −1 ; Figures 4a and 4b). The maximum (i.e., most negative) and minimum (i.e., least negative) daily ER estimates were observed in December and February, respectively. Mean ER was −469 (±135) mmol C m −2 d −1 during the snowmelt period and −239 (±92.2) mmol C m −2 d −1 during the non-snowmelt period (not shown). The standard deviation on daily ER estimates averaged 144 mmol C m −2 d −1 and had highest absolute value when the rate of ER was high but was a higher proportion of ER when it was closer to zero.
Rates of GPP were high, but large variations were observed between days and seasons. Across both years of the study, June had the highest mean GPP, while December had the lowest mean GPP (Figures 4c and 4d).  The maximum and minimum daily GPP calculations were observed in June (802 mmol C m −2 d −1 ) and December (167 mmol C m −2 d −1 ), respectively. Mean GPP was 214 (±154) mmol C m −2 d −1 during the snowmelt period and 143 (±117) mmol C m −2 d −1 during the non-snowmelt period (not shown). Throughout the study period, the standard error on the daily GPP averaged 128 mmol C m −2 d −1 . Similar to the other metabolic parameters, the absolute uncertainty on GPP was larger when GPP was larger, but values close to zero had higher relative uncertainties.
Mean NEP was negative for all months with available metabolism data (nine months), indicating net heterotrophy in Boulder Creek throughout the study (Figures 4e and 4f). The most negative monthly mean NEP was observed in December, estimated at −773 (±127) mmol C m −2 d −1 . Mean NEP during the snowmelt season was almost five times higher than mean NEP during the non-snowmelt season, with NEP during the snowmelt period averaging −239 (±161) mmol C m −2 d −1 and −49.0 (±95.4) mmol C m −2 d −1 during the nonsnowmelt period.
The metabolism model produced K 600 values that ranged from 1.33 to 43.7 days −1 and averaged 16.4 days −1 over the entire duration of the study (Table 1; Figures 4g and 4h). The mean for K 600 during snowmelt (23.5 days −1 ) was almost double of K 600 estimations in the non-snowmelt (14.3 days −1 ) seasons. Standard errors in daily K 600 estimates averaged ±8.25 days −1 . Absolute uncertainty (standard deviation) was highest when K 600 was also high, but was proportionately higher when K 600 was low, sometimes exceeding 70% of the best K 600 estimate.
Metabolism models can suffer from equifinality, a lack of uniqueness in parameter value estimates (Appling et al., 2018). The correlation between ER and K 600 was moderate for the non-snowmelt season (ρ = −0.31, p ≪ 0.01) and was higher in the snowmelt season (ρ = −0.76, p ≪ 0.01). In the complete dataset the correlation was also moderate (ρ = −0.55, p ≪ 0.01) but showed significant scatter ( Figure S4). So, despite the significant p-values, we believe the model provided interpretable results about metabolic relationships in Boulder Creek ( Figure S4).

ER Versus GPP
Examining the relationship between GPP and ER provides insight into the metabolic regime of streams and can simplify expression of the response of stream biological communities to multiple drivers (Bernhardt et al., 2018). When GPP is equal to ER (depicted by the solid black line in Figure 5a), the supply of OC

F CO2 and NEP
Most streams are net sources of CO 2 to the atmosphere, a phenomenon that can be caused by net heterotrophy or direct inputs of CO 2 from the watershed. Boulder Creek was no exception, with pCO 2 values above atmospheric equilibrium for all sensor observations and average F CO2 being positive in all months (Figure 3). Carbon dioxide flux (F CO2 ) was relatively high (97.8-1,340 mmol C m −2 d −1 ) with the highest F CO2 observed during the snowmelt period.
The relationship between F CO2 and NEP gives insight into the sources of CO 2 supporting F CO2 in Boulder Creek. A negative, one-to-one relationship (depicted by the solid black line in Figure 6a) between F CO2 and NEP suggests that net heterotrophy is the major source of excess CO 2 with little direct input of CO 2 from the watershed. Points above this line occur when the watershed contributes CO 2 to the stream which is then exchanged with the atmosphere. During the study period, a majority (99%) of F CO2 and NEP estimates were above the y = -x line, indicating that F CO2 is supported by a mixture of CO 2 produced through net heterotrophy in addition to CO 2 inputs from the watershed (e.g., terrestrial soils; Figure 6). Additionally, a strong negative correlation existed between F CO2 and NEP during snowmelt (ρ = −0.62, p < 0.001). An estimate of the percent contribution of watershed inputs of CO 2 to overall F CO2 (W CO2 ) can be made using the equation W CO2 = (100% × (1 ˗ (˗NEP/FCO 2 )). Watershed inputs were important in both snowmelt and non-snowmelt season, averaging 65% of the total CO 2 flux across the entire study. No significant difference was observed between the snowmelt (69%) and non-snowmelt seasons (63%; p = 0.59).

The Importance of Allochthonous Carbon in Boulder Creek
As observed in many other systems, day-to-day variation in stream metabolic rates was wide, with GPP ranging one order of magnitude and ER ranging nearly one order of magnitude (Table 1; Bernhardt et al., 2018;Demars, 2019;Hall et al., 2016;Hoellein et al., 2013). We observed wide variability in the daily estimated metabolic parameters, with coefficient of variation (standard deviation divided by the mean) of GPP, ER, and K 600 being 81%, 58%, and 50%, respectively, over the course of the study period. Rates of GPP in Boulder Creek were high throughout the study, with the daily average GPP in Boulder Creek, 159 mmol C m −2 d −1 , exceeding the average GPP for highly productive "summer peak" streams, 70 mmol C m −2 d −1 (Savoy et al., 2019).
We attribute these high, year-round GPP values to the hydrography of Boulder Creek at the study site. The study site exists in a hydrologically stable environment (due to water management practices) with relatively steady streamflow outside of snowmelt. This promotes the establishment of hardy primary producer communities that can metabolize year-round, even in the winter months when portions of the stream are ice-covered. Additionally, the upstream metabolic footprint of the study site and the study site itself lack a robust tree canopy cover, thereby increasing incident light to the streambed and increasing primary productivity (Bernhardt et al., 2018). This study showed a reduction in correlation between GPP and ER in the snowmelt versus non-snowmelt season ( Figures S5 and S6). Although net heterotrophy in streams is common in the literature (Battin et al., 2008;Savoy et al., 2019), some studies indicate high variability and low correlation between GPP and ER in streams (Bernot et al., 2010;Hoellein et al., 2013). For this particular study, the decoupling of GPP and ER, combined with high productivity, points to allochthonous OC fueling metabolism in Boulder Creek. Urbanization in the lower Boulder Creek watershed could be responsible for this decoupling due to consistent inputs of highly biodegradable OC (Table S1; Hosen et al., 2014). High biodegradability of allochthonous carbon would unlink GPP and ER by providing OC to heterotrophic bacteria through watershed inputs as opposed to GPP byproducts. Remineralization of urban, allochthonous carbon inputs (i.e., autochthonous carbon) at the study site is unlikely because almost all of the developed land upstream of the study site is less than 14 km upstream (calculated using USGS StreamStats; USGS, 2016) where Hall et al. (2016) estimated carbon spiraling lengths at 38-1,190 km. Additionally, total organic carbon (TOC), DOC, and carbon character does not increase substantially in the urban area upstream of the wastewater treatment plant with the exception of the snowmelt season (Murphy et al., 2003, p. 128).
ER was enhanced during the snowmelt period ( Figure 5c). The enhancement of ER during snowmelt could arise from greater mobilization of OC into the stream channel or hydrologic and biogeochemical activation of adjacent floodways and riparian areas. Although DOC and TOC concentrations are low upstream of the wastewater discharges in Boulder Creek (<3 mg C L −1 ), higher concentrations occur during snowmelt (see Chapter 5 in Murphy et al., 2003), indicating enhanced resources for heterotrophic mineralization during the snowmelt period. Additionally, DOC grab samples collected during the study period (Table S1) show, on average, higher DOC toward the beginning of snowmelt season (4.53 mg L −1 , n = 3) versus the non-snowmelt season (3.20 mg L −1 , n = 11). High flow conditions can also inundate and activate riparian zones and wetlands adjacent to streams, thereby stimulating mineralization of stored organic matter (Jacinthe, 2015), resulting in potential export of low oxygen water from these areas. This may be interpreted as enhanced stream respiration in the metabolism model used in this study.
GPP varied throughout the study and was marginally higher during the snowmelt period ( Figure 5b). The lack of reduction in GPP during increased stream discharge associated with snowmelt is similar to findings in other snowmelt-dominated streams (Ulseth et al., 2018). Increased stream discharge is believed to decrease GPP due to shading from greater water column depths and scour from sediment disturbance (O' Connor et al., 2012;Reisinger et al., 2017). However, the seasonal increases associated with snowmelt may differ, in that they can be associated with less scour or are of sufficient duration to allow re-establishment of primary producers under increased flow conditions (Ulseth et al., 2018). Spring snowmelt is also unique in that it is associated with increasing light and potentially higher nutrient inputs. These findings strengthen the concept that the ecosystem response to seasonally-driven increases in stream discharge are distinct from those caused by storms and lend support to the generality of this idea.

Seasonality of Carbon Cycling
A persistent question about stream ecosystem functioning and its relationship to F CO2 is the degree to which excess CO 2 arises from excess stream respiration (i.e., net heterotrophy) versus direct CO 2 inputs from groundwater, riparian zones, or soil. Conceptually, direct inputs of CO 2 are greatest in headwater streams, where connections to the watershed are the strongest, and decrease downstream (Hotchkiss et al., 2015;Marx et al., 2017). Our ability to pair NEP values with F CO2 for a large number of days (n = 121) enables a direct estimate of the relative contribution of NEP versus watershed inputs in sustaining F CO2 from Boulder Creek ( Figure 6). A negative correlation between NEP and F CO2 was evident during the snowmelt season ( Figure 6, ρ = −0.62, p ≪ 0.001), which supports the hypothesis that stream metabolic processes have a seasonal relationship with carbon dynamics. The moderate correlation of F CO2 to NEP during snowmelt is indicative of allochthonous OC stimulating ER, which in turn increases F CO2 . NEP and F CO2 were not significantly correlated during the non-snowmelt period (ρ = −0.12, p = 0.28); however, further knowledge about the controls on GPP, ER, and NEP in Boulder Creek, and especially their relationships to the processes controlling CO 2 flux, is needed to answer this question more fully.
Not only is the role of allochthonous OC evident in Boulder Creek through the negative relationship between NEP and F CO2 and the decoupling of GPP and ER, but there is also evidence of allochthonous CO 2 permeating the stream in the F CO2 signal. Carbon dioxide flux was higher during the snowmelt period (p ≪ 0.001), despite a lack of seasonal change in pCO 2 from the annual mean (p = 0.56; Figures 7a and 7b). K 600 was significantly higher during the snowmelt season (p = 0.003), and largely drove the increases in FCO 2 between the seasons. Short turnover time of CO 2 in Boulder Creek (on the scale of hours) throughout the year is conducive to high turnover of CO 2 , even during the snowmelt season when CO 2 storage is higher (p < 0.001; Figures 7c and 7d).

Insights Gained from F CO2 Estimates
Unusually high temporal resolution of F CO2 was achieved by combining the K 600 estimated from metabolism modeling with sensor-based pCO 2 , which provided enhanced interpretability of the biogeochemical effects of seasonality in this snowmelt-dominated mid-sized stream. Seasonal differences were particularly evident in the pairing of NEP estimates with F CO2 estimates ( Figure 7); however, the measured metabolic shift (significantly more negative NEP) during snowmelt did not translate to significant seasonal differences in NEP contribution to overall F CO2 in Boulder Creek. Despite this finding, the novel methodology of coupling modeled K 600 with high-frequency pCO 2 observations may elucidate and help to refine the view of CO 2 transport and cycling in other mid-sized watersheds.

Dynamic Aquatic-Terrestrial Linkages in the Carbon Cycle
Emerging concepts about the degree of connectivity between streams and their watersheds are often based on long-term average behavior over large spatial scales. Many have proposed that low-order, headwater streams have greater connection to the watershed resulting in more intense carbon fluxes, either as dissolved inorganic carbon, which can be emitted directly as CO 2 , or as OC, which maintains elevated ER and encourages net heterotrophy (Butman & Raymond, 2011;Butman et al., 2016;Finlay, 2011;Hotchkiss et al., 2015;Marx et al., 2017). However, the present study emphasizes the varying degrees of aquatic-terrestrial linkages that can occur at a single station depending upon seasonal flow regime. Boulder Creek acted as a low order stream during the snowmelt period and a high order stream during the non-snowmelt period with respect to CO 2 flux. Carbon dioxide flux rates increased during the snowmelt period and were likely sustained by an increase in K 600 and supplemented with increased inorganic carbon from the watershed ( Figure 6). Metabolism during the non-snowmelt period was reliant on allochthonous OC, and direct CO 2 inputs were decreased (Figures 5 and 6). Thus, the concept of the relative intensity of the processes governing stream carbon flux in headwaters versus downstream systems (Hotchkiss et al., 2015;Marx et al., 2017) is supported; however, a temporal component should also be considered, which would alter the proportion of F CO2 originating from headwater versus downstream systems within a watershed.

Conclusion
The seasonal signal shown here emphasizes the temporally dynamic nature of aquatic-terrestrial linkages, especially with respect to its influence on carbon biogeochemistry. GPP and ER were negatively correlated throughout this study, but less negative during the snowmelt period. This uncoupling of GPP and ER, along with high metabolic productivity year-round, suggests the importance of allochthonous sources of OC to maintain stream GPP and ER. Higher F CO2 during the snowmelt season further highlights the prominence and importance of allochthonous carbon inputs in Boulder Creek. This seasonal shift directly impacts the global carbon cycle by providing a direct mechanism for terrestrial OC to be metabolized to CO 2 and soilbound CO 2 to flux into the atmosphere more readily than would otherwise occur without seasonal aquatic-terrestrial linkages. The pairing of modeled metabolism estimates with F CO2 estimates was key in deciphering seasonal differences in carbon dynamics and is a novel way of estimating F CO2 on a daily timescale. The temporally dynamic nature of carbon biogeochemistry suggests that longer-term, dedicated monitoring is needed to fully resolve the connectivity and contribution of streams to the terrestrial component of watersheds and the global carbon budget. supporting this research; the Town of Erie, Colorado, for providing diversion flow data, and Hannah Podzorski and Elizabeth Whiddon for invaluable field assistance. Feedback from Rob Runkel and many anonymous reviewers helped to significantly improve this manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The authors state no conflicts of interest.