Headwater gas exchange quantified from O2 mass balances at the reach scale

Abstract Headwater streams are important in the carbon cycle and there is a need to better parametrize and quantify exchange of carbon‐relevant gases. Thus, we characterized variability in the gas exchange coefficient (k 2) and dissolved oxygen (O2) gas transfer velocity (k) in two lowland headwaters of the River Avon (UK). The traditional one‐station open‐water method was complemented by in situ quantification of riverine sources and sinks of O2 (i.e., groundwater inflow, photosynthesis, and respiration in both the water column and benthic compartment) enabling direct hourly estimates of k 2 at the reach–scale (~ 150 m) without relying on the nighttime regression method. Obtained k 2 values ranged from 0.001 h−1 to 0.600 h−1. Average daytime k 2 were a factor two higher than values at night, likely due to diel changes in water temperature and wind. Temperature contributed up to 46% of the variability in k on an hourly scale, but clustering temperature incrementally strengthened the statistical relationship. Our analysis suggested that k variability is aligned with dominant temperature trends rather than with short‐term changes. Similarly, wind correlation with k increased when clustering wind speeds in increments correspondent with dominant variations (1 m s−1). Time scale is thus an important consideration when resolving physical drivers of gas exchange. Mean estimates of k 600 from recent parametrizations proposed for upscaling, when applied to the settings of this study, were found to be in agreement with our independent O2 budget assessment (within < 10%), adding further support to the validity of upscaling efforts aiming at quantifying large‐scale riverine gas emissions.

Revisions of the global carbon cycle have indicated that rivers and other inland waters contribute substantially to the global cycling of organic carbon and emission of carbon dioxide (CO 2 ) and methane (CH 4 ) to the atmosphere (e.g., Cole et al. 2007;Battin et al. 2008;Tranvik et al. 2009;Aufdenkampe et al. 2011;Raymond et al. 2013;Hotchkiss et al. 2015;Stanley et al. 2015;Marx et al. 2017). However, due to a lack of an appropriate universal scaling for quantification of emissions, headwaters (stream order < 4) were initially not included in these assessments (see Cole et al. 2007), despite the fact they represent 17% of global riverine area (in perennial riverine systems e.g., Downing et al. 2012). Later inclusions of headwaters have proposed that small streams (order 1) might represent more than a third of the total regional emissions of CO 2 from riverine systems (Butman and Raymond 2011) and thus could be more prominent at the global scale (Raymond et al. 2013). This has highlighted the need for better constraining gas exchange at the reach scale (Trimmer et al. 2012), with the overarching goal of fine-tuning parametrizations used for large scale and global upscaling of metabolism and gas emissions.
Accurate assessments of gas exchange require quantification of the gas transfer velocity, or piston velocity (k; unit of length per unit of time −1 , e.g., m h −1 ), which physically controls the exchange of gases at the stream-air interface (see Hall et al. 2012). Thus, estimates of k are critical to the assessment of emission of greenhouse gases such as CO 2 (e.g., Hotchkiss et al. 2015;Marx et al. 2017) and CH 4 (e.g., Stanley et al. 2015;Crawford and Stanley 2016) and for the quantification of the gas exchange rate of dissolved oxygen (O 2 ) and estimates of whole-reach metabolism (e.g., Odum 1956). Stream gas exchange is most commonly assessed from total, i.e., whole-stream, metabolism estimates using open water (OW) methods based on local dissolved oxygen (O 2 ) mass balances (e.g., Odum 1956;Marzolf et al. 1994;Mulholland et al. 2001;Riley and Dodds 2013;Siders et al. 2017). The OW approach provides an integrative quantification of gas exchange R ex (in μmol L −1 unit of time −1 ) as with (O 2(sat) -O 2 ) being the O 2 saturation deficit, i.e., the difference between the O 2 concentration at saturation for the local physical conditions (O 2(sat) ) and the actual O 2 concentration in the stream water, and k 2 the O 2 gas exchange coefficient, often termed the re-aeration coefficient (k 2 or K, unit of time −1 , e.g., h −1 ). The coefficient k 2 represents the product of the depthcorrected O 2 gas transfer velocity, or piston velocity (k; unit of length per unit of time −1 , e.g., m h −1 ) and local stream areato-volume ratio, which is often approximated as the inverse of the mean water depth (in m −1 ) (see Raymond et al. 2012;Demars et al. 2015). As assessments of gas exchange critically rely on the quantification of k 2 or k, several direct and indirect methods have been developed for deriving these values. Direct measurements of these parameters are characterized by the de-gassing of a conservative tracer gas that is subsequently scaled to O 2 (e.g., Wanninkhof et al. 1990;Genereux and Hemond 1992;Reid et al. 2007;Benson et al. 2014). This approach is logistically demanding and temporal upscaling is difficult (Demars et al. 2015). Indirect measurements depend on parameterizations of the local physical characteristics of the stream reach such as depth, flow and slope (e.g., Parker and Gay 1987;Parker and DeSimone 1992 and references therein), as well as fitting of gas transfer parameters to O 2 time series based on simplified relationships (Hornberger and Kelly 1975;Chapra and DiToro 1991;McBride and Chapra 2005) or more complex modeling efforts (Holtgrieve et al. 2010;Grace et al. 2015;Appling et al. 2018). Empirical equations and parameterizations have proven to be useful as predictors of gas-exchange dynamics for specific reaches (see Genereux and Hemond 1992), but usually show large discrepancies in k 2 estimates when applied to streams with comparable hydrological characteristics (Moog and Jirka 1998;Aristegi et al. 2009;Palumbo and Brown 2014). To further constrain these parametrizations for spatial upscaling, recent studies have focused on a better characterization of the stream morphology. For example, Raymond et al. (2012) used a large collection of stream metadata to more rigorously scale k parametrizations to stream physical characteristics (i.e., flow, depth, slope, discharge).
Although k variability on short timescales, i.e., hours, has been reported (e.g., Tobias et al. 2009;, direct and indirect methods mostly focus on determining a mean value for k 2 for either the day or night, or an average for the whole day-night period, with little consideration being given to its short-term dynamics that are characteristic of most rivers. Recent application of the aquatic eddy co-variance technique (AEC) in rivers has provided robust assessment of reachscale (~150 m) benthic metabolism (Koopmans and Berg 2015;Rovelli et al. 2017) and, most recently, an "inverted" AEC approach has been used to directly quantify O 2 gas exchange in large (order 3) streams  and in an estuarine embayment (Long and Nicholson 2018). Yet, irrespective of the approach being applied to assess k, little emphasis has been given to small headwaters streams (order 1-2), despite their potential implications for the regional and global carbon cycling (e.g., Butman and Raymond 2011;Raymond et al. 2013).
The goal of our study was to provide a proof-of-concept for the derivation (values and dynamics) of k 2 and k at the reach scale in small, headwater streams (stream order 1). This was achieved by combining in situ O 2 measurements from a single station OW approach together with direct assessments of riverine metabolism using AEC, sample incubations, and groundwater inflow measurements. The proposed method was compared to the traditional nighttime regression (NR) method, and was used to characterize temporal variability in k on hourly to daily time scales. Furthermore, we evaluate k dynamics and its relation to the local physical drivers such as mean stream flow, temperature and wind patterns. The results of this study are discussed in light of parametrizations of k for regional to largescale upscaling.

Study site
This study focused on two headwaters of the Hampshire River Avon, southern England: the Chalk River Ebble (CE) and the Greensand West Avon (GA) (see Heppell et al. 2017; Supporting Information Fig. S1). For each headwater, a representa-tive~150 m reach was selected based on previous surveys of the sub-catchment morphology and geology (see "Field measurements" section). The selected CE reach (CE1 1 ; 51 1 0 41.171 00 N/1 54 0 56.309 00 W) was investigated over 3 d (25-27 April 2013) in spring, when the reach was characterized by a net outflow of stream water to the aquifer, i.e., a losing reach (see Supporting Information Fig. S2). The reach was also characterized by profuse growth (36% cover by area) of Ranunculus spp., which is widespread throughout the sub-catchment (Watson 1987). The water column of the River Ebble was characterized by low turbidity, with suspended sediment concentration limited to, on average, 37 mg L −1 (Heppell and Binley 2016a). The selected GA reach (GA2; 51 19 0 10.173 00 N/1 51 0 33.135 00 W) was investigated over 3 d from 28-30 April 2013. In contrast to the CE reach, the GA reach was characterized by a constant inflow of low-oxygenated groundwater; i.e., a gaining reach (Supporting Information Fig. S2) and only patchy macrophyte coverage (4% by area). Turbidity in the water column was 40% higher on average than at the CE reach (Heppell and Binley 2016a). The reaches are referred to as CE and GA throughout the text.

Reach scale oxygen budget
The assessment of O 2 gas exchange and stream metabolism based on a single station open-water approach (e.g., Odum 1956) relies on O 2 time series. Temporal changes in O 2 concentration, dO 2 /dt, are attributed to local whole-stream primary production (P), respiration (R) and atmospheric exchange (R ex ) as with k 2 being the gas exchange coefficient (in h −1 ) and (O 2(sat) − O 2 ) the O 2 saturation deficit. The spatial integration of a onestation approach, i.e., the integrated stretch length along the stream, is typically in the order of 1000s of m and can be quantify as 3v/k 2 , with v being the mean stream flow velocity (in m d −1 ) and k 2 (in d −1 ) (see Grace and Imberger 2006).
Rates for R and R ex are typically estimated by applying the established nighttime regression (NR) method (Hornberger and Kelly 1975). The method plots the nighttime dO 2 /dt at each time step against the O 2 saturation deficit during the night and assumes P = 0 and that respiration is light-independent. The resulting linear relationship, obtained via leastsquares regression, provides both k 2 (slope) and R (intercept). As NR-based k 2 , k NR , is only obtained at night, many studies have adopted a temperature correction term to account for temperature changes between night and day following the parameterization of Elmore and West (1961): where k 2(20 C) is the gas exchange coefficient at 20 C and k 2(T) the gas exchange coefficient at the given temperature T.
The O 2 budget (OB) approach used in the current study expands Eq. 2 by accounting for all relevant processes including (1) metabolic activity in both sediment and water, and (2) exchange with both the atmosphere and groundwater (Fig. 1). The temporal O 2 concentration was thus defined as: where F B , F WC , F K are the O 2 flux (in mmol m −2 h −1 ) associated with the benthic compartment, the water column, and atmospheric exchange, respectively. The O 2 fluxes were expressed as volumetric rates (in μmol L −1 h −1 ) by multiplying the values by the average ratio between stream area (A) and stream volume (V), which, under smooth stream water surface conditions, is approximated as the mean stream depth (d) (Demars et al. 2015).
The main strength of the proposed approach is that each flux component within the OB (Eq. 4) can be directly quantified in the field. Estimates of F B were obtained from aquatic eddy co-variance (AEC) measurements, while F WC was quantified by incubating discrete water samples in situ (see "Field measurements" section). While F B does include both streambed (F streambed ) and groundwater contributions to the O 2 budget, fluxes associated with the inflow of O 2 depleted groundwater were independently quantified as: with O 2(gw) being the groundwater O 2 concentration (in μmol L −1 ), O 2 the in-stream concentration and v gw the local groundwater inflow (in m h −1 ). The atmospheric exchange rate was quantified as: with k = k 2 (V/A) ≈ k 2 d being the O 2 gas transfer velocity, or O 2 piston velocity (in m h −1 ) and O 2(sat) the local O 2 concentration at saturation. As O 2(sat) can be quantified as a function of temperature (see Garcia and Gordon 1992), the OB approach effectively enabled a direct estimate of k 2 . To facilitate comparisons with published reference values and parameterizations from previous studies, k 2 values were also expressed as k 2(20 C) (see Eq. 3). Similarly, estimates of O 2 gas transfer velocity were also standardized to a Schmidt number of 600 2 (k 600 ) based on Jahne et al. (1987) as concentration are expressed as a function of contributions from the water column (F WC ), benthic compartment (F B ), and atmospheric exchange (F K ) following the parameterization of Eq. 4. Benthic F B includes contributions from the streambed (F streambed ) and from groundwater inflow/outflow (F GW ). Note that lateral exchange and advection processes were not considered.
where Sc T ≈ 0.00086842T 4 − 0.10115T 3 + 4.8055T 2 − 124.34 T + 1745.1 is the Schmidt number for O 2 in freshwater at the local water temperature T (see Wanninkhof 2014). Stream temperature time series were obtained from background measurements of physical parameters during the observational period at each site (Rovelli et al. 2016b).

Field measurements Site selection
The investigated reaches (CE and GA) were carefully selected to represent the sub-catchment dominant stream morphology features (e.g., stream course shape, depth, bends, riparian zone) and in-stream habitat patchiness (macrophytes, sediment type). Particular emphasis was given to ensure that each reach also fulfilled the theoretical requirements of the NR method, i.e., streams with relatively high metabolic activity and low gas exchange (see Hornberger and Kelly 1975). Field measurements were also executed under stable hydrological condition, i.e., constant flow and mean water depth over ≥ 3 d. This enabled (1) the evaluation of our local O 2 budget (OB) method and (2) comparisons with the nighttime regression (NR) method and standardized functional equations for k 2 quantification from mean hydrological parameters.

Background physicochemical measurements
Stream water column O 2 concentration, temperature, and photosynthetically active radiation (PAR) at the streambed (1 5 cm above the bottom) were monitored within each reach at 1 min intervals with a conductivity-temperature-depth logger (CTD; XR-420, RBR, Kanata, Canada), equipped with an Aanderraa O 2 optode sensor (Aanderaa, Bergen, Norway) and a 4π PAR sensor (QCP-2000; Biospherical Instruments, San Diego, U.S.). The data were also used to calibrate the O 2 sensors used for flux estimates and to define nighttime periods (PAR < 2 μmol quanta m −2 s −1 ). Wind speed and direction near the stream surface (1.5 m elevation) were recorded at 15 s intervals with a SkyWatch GEOS 11 portable weather station (JDC Electronic SA, Yverdon-les-Bains, Switzerland).

Benthic oxygen flux
The benthic oxygen flux was quantified using the AEC technique (Berg et al. 2003). Our AEC module consisted of small, lightweight, stainless steel tripod frame equipped with an acoustic Doppler velocimeter (ADV; Vector, Nortek A/S, Rud, Norway), Clark-type O 2 microelectrodes (Revsbech 1989), and submersible O 2 amperometric amplifiers (McGinnis et al. 2011). Up to two O 2 microelectrodes were positioned~0.5 cm outside the ADV sampling volume, at an inclination of~60 . Each microelectrode had a 90% response times < 0.5 s and stirring sensitivity < 0.5% (Gundersen et al. 1998). The ADV's flow measurements were obtained at 8-15 cm from the streambed with a frequency of 64 Hz and with the ADV x-axis aligned to the main flow direction. The acquired datasets were processed following the same procedure described in Rovelli et al. (2017). In short, the dataset was averaged to 8 Hz while applying data quality controls and despiking routines (e.g., Mori et al. 2007) to the dataset to remove dataset artifacts and reduce signal noise. A double coordinate rotation was applied to the flow time series to minimize the influence of horizontal motions on the vertical velocity component. The datasets were then aligned, i.e., time-shifted, to account for the relative distance between the ADV sampling volume and the microelectrode tip and for the microelectrode time constant (see Donis et al. 2015). The AEC-based turbulent oxygen fluxes (in mmol m −2 h −1 ) were estimated from timeaveraged vertical velocity fluctuations (w 0 ) and O 2 concentration fluctuations (C 0 ), as (see Berg et al. 2003) via a Reynolds' decomposition, specifically linear detrending, using the program suite Sulfide-Oxygen-Heat Flux Eddy Analysis (SOHFEA version 2.0; see McGinnis et al. 2014). The optimal detrending interval represents a trade-off between including low-frequency turbulent contributions and excluding non-turbulent contributions (McGinnis et al. 2008), and was inferred to be 5 min from cumulative averages of oxygen fluxes and friction velocity (u * ) (see Attard et al. 2014). To account for the effect of transient O 2 concentration changes in the water column between the sediment-water interface and the AEC measurement height (h), which can potentially bias F AEC (see Holtappels et al. 2013;Rheuban et al. 2014), an O 2 storage term was estimated after Rheuban et al. (2014). The benthic oxygen flux corrected for O 2 storage, F B , was defined as with dC/dt being the measured O 2 concentration gradient. The smallest area of the streambed that contributes to 90% of AEC flux, termed the "footprint area", was estimated from the sediment surface roughness parameter (z 0 ) and h following the parameterization by Berg et al. (2007). Values for z 0 were with κ being the von Karman constant (0.41), and u the average flow velocity (Wüest and Lorke 2003).

Water column activity
Water column oxygen production and consumption rates were estimated by incubating 100 mL rack-mounted glass bottles in situ over 24 h in the light and in the dark, with continuous measurements of O 2 concentration via O 2 optical fibers (4 channel FirestingO2; Pyro Science GmbH, Aachen, Germany). A set of one dark bottle, and one clear bottle were incubated at the streambed, while one bottle was incubated near the stream surface. The fourth bottle was either incubated at mid-depth at GA (depth > 0.4 m) or near the surface for the shallower CE. Three additional replicate sets were also mounted on the rack but O 2 measurements were only performed at the start and end of each incubation via Winkler titrations (Winkler 1888). Volumetric oxygen fluxes (in μmol L −1 h −1 ) were estimated from temporal O 2 concentration gradients via linear regressions for daytime and nighttime periods, respectively. No difference in fluxes were observed at the different depths, indicating a homogenously, wellirradiated water column and thus depth-independent rates. The areal oxygen fluxes F WC (in mmol m −2 h −1 ) were obtained by multiplying the values with the average stream depth. Contributions from macrophytes were deemed to be wellintegrated within benthic measurements and thus were not further addressed within the water column oxygen budget (Rovelli et al. 2017).

Groundwater influx
Dissolved oxygen flux resulting from local inflow of O 2 depleted groundwater was estimated based on measurements of the local hydraulic gradients from in-stream piezometers and porewater O 2 concentrations. Two clusters of in-stream piezometers were installed at each site,~10 m apart in the river thalweg. Each cluster comprised three piezometers screened at 20 cm, 50 cm, and 100 cm depth; the 100 cm piezometer was also fitted with narrow polytetrafluoroethylene tubing to depths of 10 cm, 20 cm, 30 cm, 50 cm, 70 cm, and 100 cm for the purposes of pore water sampling. Installation and design of the piezometers followed the description in Binley et al. (2013).
Hydraulic head was measured in all piezometers using HOBO pressure transducers (Onset Computer Corporation, Bourne, U.S.) at GA and Levelogger Edge pressure transducers (Solinst, Georgetown, Canada) at CE. Measurements were validated with manual dips on a fortnightly basis and, if needed, corrected assuming a linear drift. River stage was measured using a pressure transducer suspended within a stilling well. Falling and rising head slug tests (measurements taken by pressure transducers installed inside the piezometers) were used for computation of saturated hydraulic conductivity using the Hvorslev method (e.g., Binley et al. 2013). Vertical groundwater flux was computed using Darcy's Law, from measurements of the hydraulic gradient between the 50 cm deep piezometric head and local stream stage, combined with a weighted harmonic mean of hydraulic conductivity from slug tests carried out in the 20 cm and 50 cm deep piezometers (see Binley et al. 2013).
Pore water samples were collected from sampling tubes located on the piezometers every 2 months, from February 2014 to June 2016, using a syringe and tygon tubing. The O 2 concentration of pore water and river water was measured immediately following sample collection using a calibrated Clark-type O 2 microelectrode (50 μm tip) connected to an inline amplifier and data-logging meter (Unisense, Aarhus, Denmark) (see Heppell and Binley 2016b).

Data handling
Collected hydrological and physicochemical data were used to (1) characterize average conditions at each reach and (2) assess short term, i.e., hourly, to diel dynamics. O 2 fluxes obtained for each stream compartment were considered separately, to assess their magnitude and variability within and across the reaches, as well as within the O 2 budget context, to quantify and characterize k variability.

Reach characterization
Hydrological data and background information for each reach are summarized in Table 1. Temperature ranged from 7.8 C to 13.5 C overall, reflecting a moderate variation (< 5%) between day and night at the two sites. Daily integrated PAR (PAR 24 ) ranged between 11.9 mol quanta m −2 d −1 and 17.4 mol quanta m −2 d −1 , with peaks of up to 1630 μmol quanta m −2 s −1 . Overall, the average stream flow, based on ADV measurements (12 cm above the streambed), was in the order of 0.18-0.33 m s −1 with an associated discharge of 0.209-0.640 m 3 s −1 (Table 1) and showed no significant change between day and night, or within the observational period, indicating a stable flow. Wind speeds at CE showed irregular patterns but an overall higher wind speed during the day than at night (Fig. 2). Wind dynamics at GA, in contrast, were markedly diel, with an average wind speed of 2.4-3 m s −1 during the day and 0-0.2 m s −1 at night (Fig. 3). At both sites the water column showed well-defined diel fluctuations in O 2 concentration, with nighttime under-saturation (down to 82%) and daytime oversaturation (up to 133%), with diel O 2 concentration changes reaching~147 μmol L −1 (Figs. 2-3). The O 2 saturation deficit ranged from −108 to −96 μmol L −1 during the day, to 43-65 μmol L −1 at night. The average deficit over 24 h ranged from 2 μmol L −1 to an O 2 saturation surplus of 15 μmol L −1 .

Dissolved oxygen budget
Benthic compartment Hourly oxygen fluxes for the benthic compartment ranged between −14.0 mmol m −2 h −1 and 23.0 mmol m −2 h −1 (Figs. 2-3). At CE, oxygen fluxes followed a clear diel pattern (Fig. 2), with average daytime rates of 4.2 mmol m −2 h −1 and nighttime rates of −2.8 mmol m −2 h −1 . AEC-based benthic fluxes were estimated to cover a footprint area 26 m long and 1 m wide (~20 m 2 ), with most of the flux contribution, i.e., the region of maximum flux, located 0.5 m upstream of the AEC sampling point. At GA, a diel pattern in the benthic oxygen fluxes was also observed (Fig. 3). However, both mean fluxes at night and during the day were negative, amounting to −6.2 mmol m −2 h −1 and −3.0 mmol m −2 h −1 , respectively. The theoretical footprint area was 80 m long and 1 m wide (~70 m 2 ) with the region of maximum flux located 2 m upstream of the of the AEC sampling point.

Water column
At CE, metabolic activity in the water column was below detection (< 0.1 μmol L −1 h −1 ) and negligible for the local O 2 budget (Fig. 2). In contrast, in the more turbid, sandy GA, water   (circles), with temperature overlain (solid line). Note that the wind data were shifted by 24 h to fill the measurement gap on day 1. column activity ranged from net oxygen production rates of 2.8 μmol L −1 h −1 during the day, to oxygen consumption of −0.9 μmol L −1 h −1 at night. The average daytime rate of oxygen production was~1.2 mmol m −2 h −1 , with an average nighttime consumption rate of −0.5 mmol m −2 h −1 (Fig. 3).

Groundwater
During our study, the CE reach was losing water to the aquifer. This net local outflow of stream water had no measurable effect on in-stream O 2 concentrations (Fig. 2). At gaining reach GA, the inflow of groundwater amounted to 0.043 m d −1 on average, with a mean groundwater O 2 concentration of 63.5 μmol L −1 . Groundwater contributed to an areal O 2 concentration decrease of −0.4 mmol m −2 h −1 to −0.7 mmol m −2 h −1 (average −0.5 mmol m −2 h −1 ), representing 6-7% of the combined nighttime flux caused by metabolic activity in the benthic compartment and the water column (Fig. 3).

Atmospheric exchange
The derived rate of O 2 exchange between the stream and the atmosphere ranged from −22.3 mmol m −2 h −1 to 14.5 mmol m −2 h −1 . Estimates of k 2 based on the O 2 budget (OB) method ranged from 0.001 to 0.600 h −1, with mean values of 0.252 h −1 and 0.259 h −1 for CE and GA, respectively (Table 2). Within such range of k 2 and with mean flow velocities of 0.18-0.33 m s −1 , one-station based assessment are expected to integrate a stream length of > 7 km).
Estimates of the gas exchange coefficient k 2 using the NR method (k NR ) were 0.555 h −1 at CE and 0.434 h −1 at GA (Table 3; Supporting Information Fig. S3), and thus around a factor 2 higher than those obtained with the OB method. (Tables 2-3). At CE, we observed a significant (p < 0.01) difference between daytime and nighttime estimates of k 2 ( Fig. 2; Table 2). The average daytime k 2 (k day ) was 0.331 h −1 and about 2 times higher than k 2 at night (k night ). Similar conditions were observed at GA ( Fig. 3; Table 2), suggesting that diel changes in temperature and wind could partly explain the observed difference between day and night variability in k 2 at both sites.

Dynamics and variability of O 2 gas exchange
Temperature Gas exchange coefficient (k 2 ) and gas transfer velocity (k) are expected to be related to local stream temperature (Elmore and West 1961;Kilpatrick et al. 1989;Demars and Manson 2013). Yet, temperature has been reported to be a weak predictor of k in some settings (Tobias et al. 2009;Demars and Manson 2013), as well as in modeling studies (e.g., Correa-Gonzalez et al. 2014). For example, Tobias et al. (2009) used a modified sulfur hexafluoride (SF 6 ) tracer approach to assess variability in k (k SF6 ) on time scales of hours (in 3 h intervals) and found that k SF6 varied by 30% over a 32 h observational period and could apportion 39% of the observed variability to changes in temperature (Demars and Manson 2013).
At CE, our hourly estimates of k varied by 60%, on average, between consecutive measurements during our 3 d observational period, with temperature accounting for 46% of that variability (Supporting Information Fig. S4). As expected from previous studies (e.g., Tobias et al. 2009;Demars and Manson 2013), part of this variability could be accounted for by applying the Elmore and West (1961) temperature correction to normalize k to 20 C (k 20 ; Supporting Information Fig. S4). However, we found that temperature corrected k 20 could only explain 10% of the k-temperature relationship, and only 5% of the variability measured in k at CE. This provided further evidence that the dynamics of k could not be properly accounted for by simple functions for temperature corrections.
Instead, at CE we found that the relationship between k and temperature could be strengthened by clustering k values to fixed temperature increments, increasing stepwise from 0.5 C to 1.7 C. The best linear relationships (i.e., highest regression coefficients) were observed by clustering data in increments of 1.0-1.5 C (R 2 > 0.8), with R 2 decreasing, either side, for shorter or longer increments (Supporting Information Fig. S4) and we found the best fit to the data occurred with 1.5 C increments (Fig. 4A). As the clustering procedure partially accounts for the variability in both temperature and k, the trend in increasing fit to the data (higher R 2 ) suggested that dynamics in k and temperature, and their interaction, may occur at variable time scales. This trend was, however, not observed at GA, where temperature revealed no direct correlation with k on hourly time scales (R 2 = 0.02) or with temperature clustering (R 2 = 0.03), indicating that temperature changes in the water column did not significantly contribute to k variability at the site.
Temperature dynamics, however, may affect gas exchange without any detectable changes in bulk O 2 concentrations in the water, as recently proposed by . The authors applied the same AEC approach as in this study, but in an "upside-down" configuration, with Table 2. Estimates of k 2 (in h −1 ) using the O 2 budget (OB) method. Values for daytime, nighttime, and average over the observational period are reported as mean AE standard error, with n indicating the number of averaged data points. Standardized k 2(20 C) (in h −1 ) and k 600 (in m h −1 ) are also reported to enable better comparisons with literature studies. AEC measurements being performed near the (underside) of the stream surface (~5 cm below the atmosphere-water interface), rather than near the streambed. This enabled the quantification of O 2 flux and heat flux across the stream surface, from which they inferred local O 2 gas exchange and values for k. Their results showed that temperature fluctuations occurring just below the surface might bias the quantification of O 2 fluctuations and thus AEC-based assessments of k. The implications of the study by  suggest that other methodological approaches commonly used for the quantification of k at timescales ranging from minutes to complete diel cycles could also be impacted. The absence of any relationship between water column temperature and k at GA could thus, in part, be attributed to small-scale temperature dynamics that are not accounted for by simple temperature correction functions; although other physical parameters such as wind, for example, could also modulate k at GA.

Wind
Here the reaches were open to the wind (e.g., Supporting Information Fig. S1), with frequent wind-induced ripples being observed during the day, suggesting an interaction between flow and wind that could affect gas exchange. On time scales of hours, however, wind dynamics revealed no significant relationship with the derived variability in k (R 2 = 0.13). Similar results have also been reported from other methodological approaches, e.g., tracer release (Tobias et al. 2009) and AEC ; although both studies relied on off-site, rather than on-site wind measurements. As large variability in both k and wind might mask any general trend, wind data were temporally clustered to investigate any time lag between wind dynamics and derived k values. Using a bin interval of 1 m s −1 , which represented the dominant (> 90%) wind fluctuations at each site (Supporting Information Fig. S5), we revealed a robust positive linear relationship (R 2 > 0.7) between local wind and k 600 at each site and for both sites combined (Fig. 4B). Given the relatively low range in wind speed (0-4 m s −1 ), a linear relationship provided a suitable approximation, although for more extensive ranges in wind speed a quadratic form is more appropriate (e.g., Wanninkhof 1992;Nightingale et al. 2000;Sweeney et al. 2007;Wanninkhof 2014). Similarly to the effect of temperature, increasing or decreasing the size of the wind clustering increment resulted in a lowering of the regression coefficient, thus further supporting the argument of a temporal misalignment between the drivers of k variability and actual k dynamics.

Parametrization of gas exchange coefficients
Our study was performed under stable discharge and constant stream depth, which allowed both (1) the evaluation of our local O 2 budget (OB) method for deriving k 2 directly and then, (2) the comparison with the nighttime regression (NR) method and standardized functional equations for k 2 quantification (Table 4). The advantage of the OB method over the latter approaches is that the method is not limited to stable hydrological conditions or specific time periods. This Table 3. Estimates of k 2 using the nighttime regression method (k NR ). Values are reported with the respective regression coefficient (R 2 ) for each site. The associated plots are available in the Supporting Information (Fig. 3). The temperature-normalized values of k NR for 20 C (k NR(20 C) ), are computed based on the mean nighttime temperature at each site, while the daily NR-based k scaled for a Schmidt number of 600 (k NR(600) , in m h −1 ) was computed for the mean daily temperatures.
Night temperature ( C) contrasts with the NR method, for instance, where estimates of k 2 obtained at night might not be applicable during the day. Furthermore, the OB approach relies on in situ measurements of local O 2 mass balance variables, rather than parametrizations and scaling relationship from literature values. The method is therefore ideal for investigating Table 4. Estimates of k 600 (in m h −1 ) from established empirical equations applied to the CE and GA study sites. Input variables a : d, stream depth (m); u, flow velocity (m s −1 ); u * , friction velocity (m s −1 ); s, stream slope (m m −1 ); Q, discharge (m 3 s −1 ); and Fr ¼ u= ffiffiffiffiffi ffi gd p , Froude number with g being the gravitational acceleration constant. Note that most equations provide k 2 at 20 C (k 2(20 C) , in day −1 ; Supporting Information Table S1); k 2(20 C) values that were scaled to k 600 based on Eq. 7 and the O 2 Schmidt number for 20 C (Table 1) (Table 1). The average shear velocity, u * , computed as u * = u(C D ) 1/2 , with C D being the drag coefficient (Wüest and Lorke 2003). An average C D of 3.3 × 10 −3 was used for both sites based on further surveys of the River Avon sub-catchments (Rovelli et al. 2017). The average slope, 0.002 m m −1 , was estimated from GPS measurements during the respective field campaigns.  (Palumbo and Brown 2014). e These equations were identified to be the most suited (i.e., top performer) for the site mean depth and flow based on the suggestions of Palumbo and Brown (2014). discrepancies between the various procedures for estimating k 2 and thus strengthens the parameterization procedures applied for upscaling global estimates of outgassing coefficients for other gases such as CO 2 and CH 4 .
In this study, the NR method systematically provided among the highest estimates of k 2 for both reaches. Values of k NR were, on average, a factor of two higher than k 2 derived by the OB method (Fig. 5). The largest discrepancies, by up to a factor of three, were observed at nighttime, while, during the daytime, mean values were within reported uncertainties for this approach (Palumbo and Brown 2014). Potential issues with NR method applicability or dataset quality were deemed marginal, as both CE and GA met the requirements of the NR method (see "Methods" section), and the O 2 time series were of a high quality (Supporting Information Fig. S3). Therefore, the observed discrepancies in k 2 are attributed to a biased estimation of R derived at nighttime using the NR method (R NR ), leading to a systematic overestimation of k 2 .
Studies investigating the effect of local inflows of O 2depleted groundwater have shown that groundwater might bias the quantification of R in open water approaches and thus whole-stream O 2 budget assessments (e.g., McCutchan et al. 2002;Hall and Tank 2005;Koopmans and Berg 2015). The magnitude of such effect, which would lead to an overestimation of R, might vary according to hydrological characteristics and hydrological connectivity of each stream, ranging from negligible to substantial (see McCutchan et al. 2002). In this study, the contribution from groundwater to the O 2 budgets was quantified directly. At gaining reach GA, groundwater contributions to whole-stream estimates of R were limited to 0.5 mmol m −2 h −1 , representing < 5% of the mean R NR (15.5 mmol m −2 h −1 ); the combined R from the benthic compartment and water column activity was 6.7 mmol m −2 h −1 . At losing reach CE, contributions from groundwater were negligible but the NR method still estimated R at 8.7 mmol m −2 h −1 , which is some threefold higher than our direct measurements (2.8 mmol m −2 h −1 ).
An overestimation of R during night could also result from a temporal misalignment between actual changes in whole stream respiration and physical parameters controlling gas exchange. As shown in both datasets (Figs. 2-3), temperature decreased at night until about 6-8 a.m., and directly affected the measured saturation deficit. If k 2 is assumed to be constant, as shown from our NR regression plots (see Supporting Information Fig. S3), then the temperature-driven increase in the O 2 concentration at saturation, and associated O 2 saturation deficit, will result in an erroneous estimate of the magnitude of R NR .
It has recently been shown that k (and thus k 2 ) might display changes of up to a factor of 2-3 over timescales of 10 s of minutes to hours, even under constant hydrological conditions (e.g., Berg and Pace 2017), implying that k 2 will not remain constant over the several hours' time scale associated with the NR method. Since the NR approach considers changes in the magnitude of k 2 as a direct and instantaneous consequence of changes in R, variability in k 2 could explain the magnitude of the R NR bias obtained in this study, although the drivers of that variability, such as temperature, wind and stream discharge, need further investigation .
Hydrological parameters such as stream depth, stream flow, and stream energy (e.g., slope, shear velocity) have been used to derive k 2 (e.g.  Table 4; Fig. 5). NR-derived values were also within that range, although they ranked toward the upper end of the range. As expected from previous studies (see Aristegi et al. 2009;Palumbo and Brown 2014), estimates from each equation were highly scattered throughout the obtained range, with little consistency between the two sites. For instance, even equations that were previously identified to perform well in hydrology settings similar to our study, i.e., "top performer"  Table 4. Note that crossed equations were reported to be the "top-performer" within a specific mean depth and mean flow range from an extensive database of tracer-based k 2 values and hydrological parameters (see Palumbo and Brown 2014). Mean value and uncertainty range (dashed and dotted lines) for the hydrological parameterizations of Raymond et al. (2012) were obtained by combining mean and range from of each equation (see Table 4; Supporting Information Fig. S6). Note that to better highlight the differences between the OB, NR and common parametrizations, the equation order was re-arranged to show an incremental increase in k 600 .
equations (see Palumbo and Brown 2014), still showed large discrepancies-both within and across sites (Fig. 5). In contrast, we found that the revised functional equations provided by Raymond et al. (2012) provided well-constrained estimates, which were not only consistent within sites (see Supporting Information Fig. S6), but also between sites, with mean k 600 values (0.113 m h −1 for CE and 0.201 m h −1 for GA) within 5-7% of the local estimates from our OB method (Fig. 5).

Comments and recommendations
The aquatic eddy co-variance technique (AEC) has been shown to be an effective tool for quantifying streambed metabolism and, more recently, for the quantification of O 2 gas exchange and k from O 2 fluxes near the stream surface . A combination of the AEC technique with other traditional methods appears to be a promising next step toward better constrained assessments of gas exchange and k dynamics in headwaters, provided that careful consideration is given to (1) the site selection and (2) the representativeness of the reach within the km-sized integration of, e.g., open-water approaches. The O 2 budget (OB) approach presented in this study serves as a proof-of-concept toward this goal. While the approach is more time consuming and field demanding than traditional methods, it does enable (1) the compartmentalization of stream metabolism and (2) the assessment of short-term k variability. The study has highlighted short-term variability in k dynamics with a dampened relationship to variations of well-known physical drivers. The apparent lack of correlation could be attributed to temporal misalignments between the variability of the derived k and physical drivers, and to small-scale variability in temperature. This was clearly observed at CE, where temperature-corrected k 20 and k 600 values still revealed a marked temperature dependency, suggesting that common temperature corrections were insufficient to fully account for the observed k variability. It remains unclear whether and to what extend such small-scale variability in temperature will affect the overall gas exchange.
Although the above aspects are still poorly constrained at the local scale and on short time scales (hours), estimates of k from the most recent functional parametrizations compared well with the independent assessments of this study and should therefore be preferred over earlier parametrizations. Validations of k on local scales such as the ones presented in this study are strongly required to strengthen and add more confidence to the upscaling of k for the quantification of large-scale metabolism and global emission of climate-relevant gases such as CO 2 and CH 4 across headwaters and throughout riverine networks.