This study was conducted in four lowland headwater streams in southeastern New Hampshire and northwestern Massachusetts, United States (Figure 1). Watersheds were selected to contrast CH₄ dynamics in the two predominant land use conditions in this region (developed versus forested). These streams have been the focus of previous studies that examined nutrient and carbon cycling and the impact of land use (e.g., Wollheim et al., 2005; Wollheim et al., 2015; Wollheim et al., 2017), and were all monitored for CH₄ ebullition from May through October 2019 (Robison et al., 2021). The watersheds of two streams, Sawmill Brook (D1) and College Brook (D2), are characterized by a relatively developed, suburban landscape, while the other two, Cart Creek (F1) and Dube Brook (F2), are predominantly forest-covered. We group these streams into these two watershed land cover classes for comparison, but recognize the limitation of statistical inference with such a small sample size and mixed land use.

The watersheds range in size from 2.3 to 4.1 km², and all streams exhibit relatively shallow slopes (Table 1). Mean annual rainfall is 1,280 mm yr⁻¹ and mean annual air temperature is 8.9°C (Wollheim et al., 2017). Mean discharge ranges from 25.5 L s⁻¹ at D2 to 36.7 L s⁻¹ at D1. Stream chemistry is measured as part of the Plum Island Ecosystems LTER project (Morse and Wollheim, 2014; Wollheim et al., 2015). Stream reach characteristics were examined in a concomitant study (Robison et al., 2021). Mean reach depth at baseflow ranged from 7.2 cm at D2 to 14.5 cm at F1, and mean reach width at baseflow ranged from 1.94 m at D2 to 2.83 m at D1. All reaches exhibit shallow slopes of approximately 2 m km⁻¹ or less. Benthic substrates included sand, silt, and fine organic matter, while rocks comprised less than 10% of benthic surface area across all sites. Macrophytes were absent from all reaches in this study, eliminating the possibility of CH₄ efflux via plant-mediated transport.

Water samples for dissolved CH₄ analysis were collected at each stream using 60 ml syringes fitted with three-way stopcocks. Samples were collected at variable frequencies ranging from multiple samples a week to roughly weekly between June 1 and October 31, 2019. On average, a sample was collected every 6 days at each site. Syringes were rinsed with stream water prior to sample collection. To collect water samples, syringes were filled with approximately 60 ml of stream water from 5 to 10 cm depth below the stream water surface. Syringes were cleared of air bubbles by inverting and expelling bubbles and water until 30 ml of sample water remained. Samples were stored on ice until returned to the laboratory within 6 h. In the laboratory, 30 ml of ambient air was added to each syringe to achieve a 1:1 ratio of sample water to air. Syringes were then shaken for 2 min to equilibrate gases between water and headspace (Magen et al., 2014). The water was then dispelled from the syringe, and the remaining headspace gas was saved for analysis. If the gas samples were not analyzed immediately, they were stored in evacuated glass vials sealed with a rubber septum until analyzed.

Samples were analyzed for CH₄ concentration in the Trace Gas Biogeochemistry Laboratory at the University of New Hampshire. The CH₄ concentration in parts per million by volume (ppmv) was determined using a Shimadzu Gas Chromatograph Flame Ionization Detector. Concentration was standardized using the average area response of ten injections of a standard CH₄ mixture (Northeast Airgas, 2.006 ppmv or Maine Oxy, 1,000 ppmv) to determine instrument precision (Frolking and Crill, 1994). If multiple samples were collected from a single site, the mean measured concentration was used.

Diffusive fluxes of CH₄ to the atmosphere (F_CH4, mmol CH₄ m⁻² d⁻¹) were calculated as: F CH 4 = k CH 4 × Δ CH 4 , (1) where k_CH4 (m d⁻¹) is the gas transfer velocity for CH₄ and ΔCH₄ is the difference in CH₄ concentration (g m⁻³) between the water and the air corrected for Henry’s Law. We used a constant air CH₄ concentration of 1.94 ppmv based on measurements at the nearby Global Monitoring Laboratory site at the Isle of Shoals, New Hampshire during the monitoring period in 2019 (Earth System Research Laboratory, G. M. D., 2021). The gas exchange rate for CH₄, k_CH4, was calculated as: k CH 4   = ( SC CH 4 600 ) 1 / 2 /  k 600 , (2) where SC CH 4 is the Schmidt number for CH₄ at a given water temperature and  k 600 (m d⁻¹) is the gas transfer velocity standardized for a Schmidt number of 600. Without direct measurements of gas transfer velocity across a wide range of flow conditions, we estimated  k 600 for each reach based on the relationship following Raymond et al. (2012): k 600   = 1162   ± 192 × S 0.77   ± 0.028 × V 0.85 ± 0.045 , (3) where V is the water velocity at time of sampling (m s⁻¹) and S is the channel slope (unitless). The mean slope of each stream reach was found using StreamStats (Ries et al., 2017). Mean daily discharge (Q; L s⁻¹) was available from the long-term monitoring projects (Morse and Wollheim, 2014; Wollheim et al., 2017). Predictive relationships for V from Q were made using the equations for scaling stream geometry at a site (Knighton, 1998): V  = 0.287 × Q 0.4 . (4)

Two or fewer direct measurements of the gas exchange velocity were also made at each site, although generally under relatively low flow conditions (< 50 L s⁻¹). These measurements were made using argon as a conservative gas tracer (Hall and Ulseth, 2020). We used these measurements to roughly evaluate the accuracy of our estimates of k₆₀₀ (Supplementary Figure S1).

Ebullitive CH₄ fluxes were measured during the same time frame as diffusive flux by Robison et al. (2021). Briefly, stationary bubble traps were deployed in triplicate at three to four locations in each stream. Traps were visited at least weekly from June 1 to October 31, 2019, during which the volume of gas collected in each trap was measured and collected for analysis of CH₄ concentration. The flux at each trap was then calculated as the mass of CH₄ emitted over the observation period normalized per unit area under the trap (mmol CH₄ m⁻² d⁻¹). For this study, mean rates of ebullitive CH₄ flux per measurement period at each stream were used, based on the total number of funnels at a stream.

Gas samples were collected for CH₄ isotopic analysis in the first week of August 2019. Additional samples from F1 and D1 from the first week of September 2018 were also included in analysis. Dissolved CH₄ samples from the surface water were collected in the same manner described above. Additional benthic gas samples were collected by physically disturbing the sediment and collecting released bubbles. A minimum of three sediment gas samples and three surface water dissolved gas samples were collected at each site. All samples were stored in evacuated glass vials sealed with a rubber septum for analysis.

Samples were analyzed for δ¹³C-CH₄ and δD-CH₄ using an Aerodyne dual tunable infrared laser direct absorption spectrometer (TILDAS; Aerodyne Research Inc, Massachusetts, United States ) at the University of New Hampshire. These instruments use high resolution infrared spectrometry to quantify trace gases such as CH₄ (Mcmanus et al., 2011; Nelson and Roscioli, 2015). The TILDAS used here is configured with two 8 µm quantum cascade lasers (QCLs, Alpes Lasers, Switzerland) and a 200 m multipass absorption cell to simultaneously monitor ¹²CH₄, ¹³CH₄, and CH₃D. The instrument was regularly calibrated with three standard tanks with known isotopic mixing ratios for ¹³CH₄/¹²CH₄ and CH₄/CH₃D. The isotopic composition of the standards was determined using an Aerodyne calibration system. The spectroscopic isotope ratios of four Isometric (now Airgas) CH₄ standards were measured at diluted concentrations ranging from < 1 to 12 ppmv. Keeling plot analysis was used to determine the relationship between the spectroscopic and standard isotope ratios of the Isometric standards. The linear relationship was made to the corresponding measured spectroscopic isotopic ratios of the UNH calibration tanks. The TILDAS is configured with an automated sampling system designed to measure small (≤ 5 ml) injections of high concentration samples by diluting them within the instrument to a target CH₄ mixing ratio of eight ppmv with ultra-zero air. The instrument precision was 0.1‰ for δ¹³C-CH₄ and 3‰ for δD-CH₃D at the target CH₄ mixing ratio. Samples with very low concentration (<0.01 mM CH₄) can result in erroneous isotopic measurement, and these datapoints (five of 59 samples) were removed prior to statistical analysis.

The percent of CH₄ oxidized between sediment generation and surface water evasion was estimated using the equation: % oxidized =   δ 13 C Sample − δ 13 C Source ( α − 1 ) × 1000 × 100 , (5) where δ 13 C S a m p l e is the δ¹³C-CH₄ value found in the surface water sample, δ 13 C S o u r c e is the δ¹³C-CH₄ value found in the sediment gas sample, and α is the isotope fractionation factor of CH₄ oxidation. We used two α values to include some uncertainty of isotopic fraction in these stream environments, 1.033 and 1.025 (Tyler et al., 1997; Zhang et al., 2016). This follows previous work estimating the efficiency of CH₄ oxidation in the Amazon River basin (Sawakuchi et al., 2016).

Sediment cores were collected from the streams for 16S rRNA analysis. A modified Multi Stage Soil Core Sampler (AMS Inc, Idaho, United States ), consisting of a stainless steel cylinder and a plastic liner (5 cm diameter), was manually driven into the stream sediments using a sliding weight stand (Wik et al., 2018). Sediment cores were collected the first 2 weeks of July 2019 and processed in the field. Cores were collected at patches in the stream near each set of bubble traps as described in Robison et al. (2021). Four patches at each stream were sampled, except D2 where only three patches were sampled. At each patch, triplicate cores were collected. Sediment collected in the top 2 cm from each core was combined into a composite “surface” sediment sample in a Whirl-Pak sample bag. Similarly, sediment collected between 9 and 11 cm depth was combined into a composite “subsurface” sample. As a result, a total of four composite surface and four composite subsurface samples were collected at each stream except D2, where only three composite samples of each were collected. All tools used to collect cores, including the coring equipment, were sterilized with 70% ethanol between samples and sediment layers to minimize contamination. The depth of 9–11 cm was chosen as the subsurface sample depth because buried rocks or clay layers prevented collection of sediment at depths greater than 12–15 cm in some places. Thus, while 9–11 cm does not represent the maximum or mean depth of the sediment in portions of these streams, it represents a depth layer that could be collected across all sites and is deep enough where complete depletion of oxygen is likely (Crawford et al., 2014). Sediment samples were frozen on-site with liquid nitrogen, stored on ice, and returned to the laboratory within 4 h. Samples were stored in a −80°C freezer until analyzed.

Genomic DNA was extracted from sediment samples using a Qiagen DNeasy PowerSoil kit (Qiagen, Hilden, Germany) with minor changes to the manufacturer’s protocol (see Doherty et al., 2020). DNA was amplified by polymerase chain reaction (PCR) using the primers 515f–806r of the V4 region of the 16S rRNA gene to profile bacterial and archaeal communities (Apprill et al., 2015; Parada et al., 2016). Each reaction contained 6 µl DreamTaq Hot Start Green (Thermo Fisher Scientific, Waltham, MA, United States ), 2.6 µl sterile water, 0.7 µl forward primer (5 µM), 0.7 µl reverse primer (5 µM), and 2 µl template DNA (10× diluted). DNA amplification was performed using a T100 Thermal Cycler (Bio-Rad, Hercules, CA, United States). The 16S rRNA conditions were: enzyme activation at 95°C for 3 min; followed by 35 cycles of denaturation at 95 °C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 60 s; then a final extension at 72°C for 12 min.

Gel electrophoresis was used to confirm the presence of the PCR product, which was then quantified using a Quant-iT dsDNA High-Sensitivity Assay Kit and a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Massachusetts, United States). Resulting concentrations ranged from 0.5 to 28 ng μL⁻¹, with sandier sediment samples typically having lower concentrations. PCR products were sent to the Hubbard Center for Genome Studies (University of New Hampshire, New Hampshire, United States) for sequencing by Illumina HiSeq2500 (250 bp reads) with Rapid Run^© SBS V2 chemistries (Illumina, San Diego, CA, United States). Sequence reads were demultiplexed using CASAVA.

Sequences were analyzed using QIIME 2 (version 2020.1; Bolyen et al., 2019) on the Premise high performance computing cluster at the University of New Hampshire. Cutadapt was used to remove primers (Martin, 2011). Using DADA2 (Callahan et al., 2016), sequences were truncated at the length where the median Phred score (quality score) fell below 30 (200 bp for forward read; 225 bp for reverse read) and unique sequences were designated as amplicon sequence variants (ASVs). Each sample was then rarefied to 3,400 sequences per sample. Rarefication depths were chosen to ensure at least three samples remained for each depth at each stream. Taxonomy was assigned to ASVs using the SILVA database (release SILVA 138 SSU; Quast et al., 2013; Glöckner et al., 2017). The relative abundances of taxa were calculated after removing chloroplasts, mitochondria, and taxa present in less than 5% of samples.

In this study, we focused primarily on taxa associated with CH₄ production and CH₄ oxidation. Following Hough et al. (2020), we identified methanogens based on classification into the archaeal orders: Methanocellales, Methanobacteriales, and Methanomicrobiales, and the families Methanosarcinaceae and Methanosaetaceae, of order Methanosarcinales (Evans et al., 2019). Similarly, we identified methanotrophs as those from the bacterial orders Methylococcales and Methylomirabilales and the archaeal family Methanoperedenaceae, order Methanosarcinales (Smith and Wrighton, 2019). Taxa from the order Methylococcales are classical aerobic methanotrophs, while those from the order Methylomirabilales are capable of CH₄ oxidation following reduction of nitrate or nitrite (Wu et al., 2011); that is, they are capable of generating oxygen for CH₄ oxidation when in an anoxic environment. Finally, Methanoperedenaceae are a methanotrophic archaea family (formerly known as ANME-2d) known to be capable of CH₄ oxidation via alternative electron donors like nitrate, iron, or manganese (Ettwig et al., 2016; Leu et al., 2020). Despite the indefinite nature of assigning organismal function using 16S sequencing, traits associated with CH₄ cycling are typically conserved and well defined (Martiny et al., 2015), allowing for confidence in these assignments. The relative abundance of individual samples is calculated, from which a mean and standard deviation from all samples collected at a sediment depth in a stream is determined.

All calculations and statistical analyses were performed in MATLAB and Statistics Toolbox Release 2020a (The MathWorks, Inc, Massachusetts, United States), and the significant level was set at α = 0.05. Uncertainty in the dissolved CH₄ concentration for each site was determined by the ratio of the standard error (SE) of the measured concentrations to the median concentration (Hojo and Pearson, 1931): % uncertainty  = 1.253 ×   Std Err median CH 4 × 100 , (6) which resulted in a unique percent uncertainty at each site. Uncertainty in the gas exchange rate was determined by the standard deviation included in the factors of Eq. 3. Simple linear interpolation was used to estimate the rate of diffusive efflux on days without a dissolved CH₄ concentration measurement. Uncertainty in the diffusive CH₄ flux was then estimated by propagating the uncertainty from dissolved CH₄ concentration and gas exchange rates. We summarized CH₄ concentration, gas exchange, and rates of diffusion as medians because of the non-normal distribution of the data. Differences in CH₄ concentration, gas exchange, and diffusive emissions between streams were analyzed using the nonparametric Wilcoxon rank sum test. Trends in CH₄ concentration and diffusion across time were analyzed using the nonparametric Kendall rank correlation coefficient. Isotopic and microbial data were summarized as means because of their small sample sizes and unknown distributions. Differences between streams, depths for microbial samples, and type of sample for isotopes (i.e., benthic or surface water) were compared using two-sample t-test or one-way ANOVA.

Measured concentrations of dissolved CH₄ in surface water were above saturation (roughly 0.003 µmol CH₄ L⁻¹ at 10°C) in all collected samples, with an overall median ±SE of 1.1 ± 0.4 µmol CH₄ L⁻¹ (Table 2; Figure 2A). F2 and D2 had the highest median concentration of CH₄ (1.6 ± 0.9 µmol CH₄ L⁻¹) and F1 had the lowest (0.9 ± 0.3 µmol CH₄ L⁻¹). There was no statistical difference in the median dissolved concentration of CH₄ between watershed land cover classes (p > 0.05). We also observed no clear pattern in the dissolved concentration across time at any individual stream (p = 0.54). The median k_CH4 across all sites was 4.8 ± 0.4 m d⁻¹, ranging from 3.4 ± 0.3 m d⁻¹ at F2 to 5.9 ± 0.6 m d⁻¹ at D1 (Table 2; Figure 2B). The watershed land cover classes did differ in the estimated gas exchange rates (p < 0.01), where the two streams draining developed landscapes, D1 and D2, exhibited significantly higher k_CH4. Calculated rates of gas exchange were similar to limited direct measurements in each of the streams (Supplementary Figure S1) using short-term continuous injections of a volatile gas tracer (Hall and Ulseth, 2020). Median diffusive CH₄ flux was 5.7 ± 2.3 mmol CH₄ m⁻² d⁻¹ overall (Table 2; Figure 2C). In individual streams, the median diffusive flux ranged from 4.1 ± 2.1 to 7.8 ± 4.8 mmol CH₄ m⁻² d⁻¹ at F1 and D2, respectively, with the developed watershed class exhibiting a higher rate compared to the forested watershed class (p < 0.01). We did not observe any clear temporal pattern in diffusive CH₄ emissions across our monitoring period at any individual site (p > 0.05). High diffusive CH₄ fluxes were more strongly associated with high measurements of dissolved CH₄ than with high gas exchange rates (Figure 2). The maximum daily rate of diffusive CH₄ efflux was 70.8 ± 31.4 mmol CH₄ m⁻² d⁻¹ at D2.

Diffusive fluxes of CH₄ were compared to ebullitive fluxes of CH₄ from Robison et al. (2021) during the same period of record (Table 2; Figure 2D). Across all four streams, the median diffusive CH₄ flux (5.7 ± 2.3 mmol CH₄ m⁻² d⁻¹) was greater than the median ebullitive CH₄ flux (0.6 ± 0.1 mmol CH₄ m⁻² d⁻¹). Diffusive emissions of CH₄ comprised approximately 90% of total CH₄ emissions from these four streams over the entire period, ranging from 78% at F2 to nearly 100% at D2. Ebullitive CH₄ flux exceeded diffusive flux only at D1 for short periods (Figure 2D). Diffusion accounted for the majority of CH₄ emissions at the other three sites during the entire monitoring period.

Benthic gas and surface water gas samples exhibited clear differences in δ¹³C-CH₄ and δD-CH₄ values (Table 2; Figure 3). Generally, there was more variability between the type of sample (i.e., benthic or surface water) than between streams. Benthic gas δ¹³C-CH₄ and δD-CH₄ values averaged across all streams were −65.9 ± 2.0‰ and −300 ± 18‰, respectively. There was some variation in the benthic δ¹³C-CH₄ values between streams, where D1 exhibited a significantly higher mean δ¹³C-CH₄ (−63.1 ± 0.6‰) than F1 (−67.8 ± 0.5‰, p < 0.01) and F2 (−66.5 ± 0.4‰, p = 0.02). More variation was noted in measured δD-CH₄ values, with F2 having the highest mean value of −257 ± 16‰ and D2 having the lowest value of −318 ± 16‰, with only these two sites being significantly different (p = 0.01). Dissolved CH₄ δ¹³C-CH₄ and δD-CH₄ values averaged across all four streams were −52 ± 1.5‰ and −261 ± 29‰, respectively. These were slightly more varied across streams than the benthic samples. For example, δ¹³C-CH₄ values ranged from −49.5 ± 0.6‰ at D1 to -55.0 ± 2.6‰ at D2 in dissolved CH₄ samples, although no site was significantly different from the others. Meanwhile, dissolved δD-CH₄ ranged from −257 ± 26‰ at F2 to −290 ± 42‰ at D1, with only these two sites being significantly different (p = 0.04). There was no difference in measured isotopic values between watershed land cover classes when comparing either δ¹³C-CH₄ or δD-CH₄ between benthic gas and dissolved samples.

Dissolved CH₄ in surface water was more ¹³C- and D-enriched relative to the sediment gas overall (p < 0.01; Figure 3). This amounted to a 14 ± 1.5‰ and 39 ± 29‰ shift towards the heavier isotope for δ¹³C-CH₄ and δD-CH₄, respectively. The δ¹³C-CH₄ values for dissolved CH₄ were all significantly higher than sediment CH₄ at each stream, and the δD-CH₄ values were significantly higher than sediment CH₄ at F1 (p = 0.01) and D1 (p = 0.04). Given this enrichment in δ¹³C-CH₄ values, we estimated the mean proportion of dissolved CH₄ oxidized was 42.1 ± 2.1% when using α = 1.033, and 55.6 ± 2.8% when using α = 1.025 (Table 2). The percentage of dissolved CH₄ oxidized in individual streams ranged from a minimum of 33.3 ± 3.1% at D2 using α = 1.025 to 65.2 ± 2.6% at D1 given α = 1.033.

Taxa identified as methanogens and methanotrophs were detected across all four streams (Figure 4A, Supplementary Table S1). Methanogen taxa were detected in 10 of 15 surface sediment samples, and 14 of 15 subsurface samples (Figure 4B). These included the archaeal orders Methanosarcinales (families Methanosaetaceae and Methanosarcinaceae), Methanobacteriales (family Methanobacteriacea), Methanomicrobiales (families Methanomicrobiaceae, Methanoregulaceaea, and rice cluster II), and Methanocellales (family Methanocellaceae). As a percent of the total community, methanogens comprised roughly 0.1% of the surface sediment community and 0.6% of the subsurface community, on average. Methanogens comprised a larger percentage of the microbial community in the forested watershed streams than in the developed streams (p = 0.04).

Taxa classified as methanotrophs were detected in 14 of 15 surface samples and all subsurface samples (Figure 4C). Taxa from the bacterial orders Methylococcales (families Methylococcaceae and Methylomonaceae), Methylomirabilales (family Methylomirabilaceae), as well as the archaeal order Methanosarcinales (family Methanoperedenaceae) were detected. Methanotrophs comprised a larger fraction of the total microbial community in all samples compared to methanogens, but this may be a result in part of the bias of our primers for bacteria over archaea (Walters et al., 2016). As a percent of the total community, methanotrophs comprised on average 4.3 and 3.6% in surface and subsurface samples, respectively. Aerobic methanotrophs, i.e., Methylococcales, were significantly more abundant in surface sediment samples (p = 0.02), while taxa capable of CH₄ oxidation in anoxic conditions (Methylomirabiales and Methanoperedenaceae) were higher in relative abundance in subsurface samples than in surface samples (p = 0.02). Additionally, methanotroph relative abundance was greater in the two streams draining developed watersheds than the two draining forested watersheds (p < 0.01).

Diffusive emissions were the dominant efflux pathway of CH₄ in these four lowland stream reaches. The median diffusive flux is comparable to a summary of 385 streams and rivers (8.22 ± 25.50 mmol CH₄ m⁻² d⁻¹; Stanley et al., 2016). Total CH₄ emissions were greater than the median of a recent analysis of streams and rivers globally (0.9 mmol CH₄ m⁻² d⁻¹), but within the 90% range (0.0–35.8 mmol CH₄ m⁻² d⁻¹; Rosentreter et al., 2021). When compared to lentic waterbodies, the rate of emissions from these four streams is similar on a per area basis (e.g., Holgerson and Raymond, 2016; Sanches et al., 2019). For example, in a variety of lakes and ponds of northern latitudes, total CH₄ emissions ranged from 2.8 to 12.5 mmol CH₄ m⁻² d⁻¹ (Wik et al., 2016b).

The efflux of CH₄ from small streams may be a significant, yet overlooked pathway of carbon fluxes from inland waters. The total aquatic carbon flux for inland waters across North America is estimated to be 24 g C m⁻² yr⁻¹ (Butman et al., 2018), where the flux is normalized for total land area. However, this flux does not include CH₄, as scaling was limited by significant data gaps at the time of estimation. Following scaling methods published previously (Robison et al., 2021), the streams in this study emit roughly 0.08 g CH₄-C m² watershed area yr⁻¹ when both diffusion and ebullition are considered, less than 1% of the total estimated for inland waters. If the greater radiative forcing of CH₄ relative to CO₂ is taken into account (Saunois et al., 2020), the CO₂-equivalent rate increases to 2.2 g C m² watershed area yr⁻¹. Thus, ignoring small stream methane emissions may underestimate watershed carbon emissions by approximately 10%, a significant factor.

Relatively few studies exist with which to compare diffusion and ebullition in lotic systems (Supplementary Table S2). Summarizing 26 streams and rivers which considered both pathways, diffusion delivered 68% of total CH₄ emissions on average, but this ranged from 14 to 94% (Stanley et al., 2016). However, it appears smaller streams tend to favor diffusive emissions. For example, diffusion comprised roughly 90% of CH₄ emissions in a small stream in Wisconsin, United States (Crawford et al., 2014). In contrast, studies of larger rivers in the Amazon (Sawakuchi et al., 2014), on the East Qinghai–Tibet Plateau (Zhang et al., 2020), and in urban China (Wang et al., 2021) found diffusion accounted for 64, 21, and 20% of total CH₄ emissions, respectively. Diffusive fluxes generally account for a smaller fraction of the total CH₄ budget in lentic systems too, e.g., 48% for northern lakes and ponds (Wik et al., 2016b).

Land use may affect diffusive CH₄ emissions, where development promotes higher rates of emission. This pattern follows previous studies which linked low oxygen and high organic carbon delivery to higher rates of CH₄ production in urban streams (Stanley et al., 2016; Wang et al., 2021). With only four streams studied and mixed land use within each watershed, the strength of this comparison is narrow. Individual samples from each stream exhibited relatively high concentrations of dissolved CH₄. The median and maximum of dissolved CH₄ concentrations were similar to those found in a yearlong survey of gas concentrations in streams in New Hampshire (0.63 and 13.5 µmol CH₄ L⁻¹; Herreid et al., 2020). Our dataset precludes analysis of what causes these high concentrations, but we have no reason to discard them from our analysis as they fall within ranges observed in streams and rivers locally and globally (Stanley et al., 2016).

The estimation of diffusive CH₄ efflux is limited by the relatively sparse measurements of dissolved CH₄ concentrations and the estimation of gas exchange. The similarity of direct gas exchange measurements to estimates used in this study suggests this latter uncertainty is modest during low flow conditions (Supplementary Figure S1). Similarly, while our results likely suffer in accuracy at fine temporal scales (i.e., daily), it does not necessarily indicate the summarized results are inaccurate. Uncertainty due to sampling frequency on estimating seasonal rates of CH₄ emissions in lakes has been examined previously (Wik et al., 2016a); here, a minimum of 17 measurements of dissolved CH₄ concentration were needed to accurately represent seasonal emission patterns from lakes. The minimum number of measurements of dissolved CH₄ concentration in this study is 25. While the standard for lakes may not be the same for lotic ecosystems, we believe our results reflect the typical seasonal behavior at these streams in which diffusion is the dominant emission pathway for CH₄.

The isotopic and microbial results provide evidence for active CH₄ processing in these streams, with production and oxidation ubiquitous. For example, the presence of ebullition and the detection of methanogens suggests active CH₄ production in the sediment of all four streams. The δ¹³C-CH₄ and δD-CH₄ values of the sediment-derived gas suggest a combination of acetoclastic (−60‰ to −50‰ δ¹³C-CH₄ and −400‰ to −250‰ δD-CH₄) and hydrogenotrophic (−110‰ to −60‰ δ¹³C-CH₄ and −250‰ to −170‰ δD-CH₄) methanogenesis (Whiticar, 1999; Chanton et al., 2005), which is common in aquatic environments (Shelley et al., 2015; Wik et al., 2020). The presence of methanotrophs and an enrichment in both δ¹³C-CH₄ and δD-CH₄ values between benthic and surface samples suggests CH₄ oxidation was prevalent across the four streams as well, representing a potentially significant control on emissions. The estimated oxidation proportions are similar, but generally lower than the 57–82% range estimated for large rivers in the Amazon River basin (Sawakuchi et al., 2016) or the 57–100% estimated for lakes in Sweden (Bastviken et al., 2002). The lower oxidation proportion in the streams could be a result of relatively faster turnover of the water column. That is, while the efficiency of CH₄ oxidation could be similar across all aquatic ecosystems, relative fast mixing in small streams may limit the time for CH₄ oxidation to occur. Regarding microbial evidence for CH₄ production and oxidation, detection of microbial taxa via 16S rRNA sequencing can recover dormant or dead microbes (Oliver, 2005), thus the detections do not directly indicate activity. However, combined with the magnitude of CH₄ fluxes and isotopic values, we can confidently infer the presence of microbial CH₄ production and oxidation in the streams.

Alternative explanations for the change in isotopic values are possible. The difference could result from dissolved CH₄ inputs from a different source in the watershed, e.g., the riparian zone (Crawford et al., 2013). If the balance of isotopic values of this other source is distinct from that in the stream sediments, the isotopic values of the dissolved CH₄ in the streams may be confounded by the mixed sources. Additionally, diffusion exerts a fractionation effect (< 3‰ δ¹³C-CH₄) on CH₄, where lighter CH₄ is preferentially emitted (Knox et al., 1992). This fractionation is much smaller than that of oxidation and the time for diffusion is relative fast in these streams, limiting the impact of this effect. Because isotopic and microbial sampling occurred in August and July, respectively, the results of our oxidation analysis may reflect the conditions of this warmer period of the year. For example, CH₄ isotopes in aquatic ecosystems may be heaviest in the summer (Atkins et al., 2017), possibly as a result of more active microbial oxidation (Shelley et al., 2015). If so, this would indicate we may have measured conditions under which oxidation is most efficient. However, considering the isotopic results in the context of the detection of methane oxidizing microbes, we are confident CH₄ oxidation is prevalent in these streams with the greatest uncertainty around the temporal dynamics of oxidation efficiency.

The dominance of diffusive CH₄ emissions and the relatively lower proportion of CH₄ oxidation appear to distinguish small streams from other aquatic ecosystems. Identifying the factors that drive this difference is thus key to understanding how CH₄ dynamics may differ in small streams. We propose two characteristics that are fundamental in differentiating small streams from larger rivers and from lentic systems: flow and depth. We hypothesize limited bubble formation and CH₄ oxidation results from the unique physical and biogeochemical conditions of small streams caused by increased water exchange between the water column and sediments (hyporheic flow) and short residence time of CH₄ in the water column (gas exchange).

Flow provides lotic ecosystems with consistent turbulent mixing, and the relatively shallow nature of small streams propagates this turbulence to the sediment–water interface. These likely limit bubble formation broadly via changes to sediment quality and depth. Because bubble formation is greatly affected by the grain size and physical structure of aquatic sediments (Liu et al., 2016), processes like turbulence which affect the sediment will in turn affect ebullition. For example, the erosional nature of streams (Knighton, 1998) may inhibit ebullition by increasing mean grain size or reducing deposition of organic matter, both of which have been shown to limit ebullition in streams (Crawford et al., 2014; Bodmer et al., 2020). Turbulence also enhances exchange at both the sediment–water and air–water interface within streams (Boano et al., 2014), affecting both bubble formation and oxidation rates. Advective water inputs into and through stream sediments via hyporheic flow (Packman and Salehin, 2003) create patches of oxic and anoxic conditions in close proximity (MahmoodPoor Dehkordy et al., 2019; Nelson et al., 2019). The CH₄ produced in anoxic zones can be transported through sediments by this advective flow (Sobczak and Findlay, 2002), limiting bubble formation further, exposing CH₄ to oxidizing conditions, and decreasing the residence time of CH₄ within the sediment. This advective exchange may ultimately limit the overall potential for CH₄ oxidation. For example, in Amazonian rivers, the extent of CH₄ oxidation was positively related to residence time in surface sediments (Sawakuchi et al., 2016). Oxidation within the water column is also limited by higher rates of gas exchange in streams and rivers compared to lentic ecosystems (Looman et al., 2021), as does the shallower depth of streams compared to larger rivers (Raymond et al., 2012). For example, 51–100% of the CH₄ in lakes originating in sediments in deeper portions of lakes was oxidized in the water column, while only 24–40% of the CH₄ from shallower areas was oxidized due to the relatively shorter transit time (Bastviken et al., 2008).

We propose a conceptual model of CH₄ production, transport, oxidation, and emission unique to small streams, where CH₄ is emitted primarily through diffusive fluxes and oxidation is relatively limited (Figure 5). In this model, 1) CH₄ is primarily produced by methanogens in anoxic sediments. While bubble formation is generally limited in small streams, 2) some CH₄ is emitted via ebullition. The bulk of CH₄ production, plus inputs from the upslope catchment via groundwater flow, remains dissolved. In small streams, dissolved CH₄ 3) diffuses through the sediments entering hyporheic flow paths, 4) although small amounts of vertical diffusion through the sediments may also occur. Hyporheic flow, which remains primarily anoxic in many systems, moves dissolved CH₄ back to the water column. We propose 5) oxidation primarily occurs near the sediment–water interface, where the redox gradient is greatest and the highest abundance of methanotrophs are detected. In the surface water, CH₄ oxidation is likely limited due to relatively high gas exchange rates. After oxidation, 6) the remaining dissolved CH₄ emitted via diffusion is enriched isotopically relative to that emitted via ebullition. The relative importance of each of these processes will depend on stream characteristics such as slope, water velocity, water depth, and sediment characteristics (i.e., depth, grain size) that influence the degree of hyporheic exchange, hyporheic residence time, redox conditions, and microbial activity. Some of the proposed mechanisms described in this conceptual model were not directly measured (e.g., hyporheic exchange) and require further exploration and testing. We specifically urge the examination of CH₄ cycling and transport within the benthic environment in relation to subsurface flow and redox conditions.

One important implication of this conceptual model is how streams are considered in global CH₄ models (Dean et al., 2018), particularly with regards to the isotopic values of emitted CH₄. For example, given both the imbalance of emission pathways in this study, as well as the isotopic values of ebullitive and diffusive gases, we calculate mean isotopic values of −55.7 ± 2.1‰ δ¹³C-CH₄ and −260 ± 36‰ δD-CH₄ for the total CH₄ emitted. The δ¹³C-CH₄ reflects relatively heavy values compared to global averages of approximately −62‰ δ¹³C-CH₄ used for aquatic systems (Schwietzke et al., 2016). Considering streams and rivers account for roughly 8% of the global aquatic CH₄ budget (Rosentreter et al., 2021), this difference could be enough to alter the global isotopic signature used for aquatic CH₄ sources.

Very few studies which combine measurements of CH₄ isotopes and emission pathways are available in lotic systems, limiting comparison of our results (Supplementary Table S2). The most extensive data in rivers comes from the Amazon, where consideration of both ebullitive and diffusive emissions provides mean δ¹³C-CH₄ values for total emitted CH₄ of −49.6‰ (Sawakuchi et al., 2016) and −38.4‰ (Sawakuchi et al., 2021). In these large tropical rivers, diffusion comprises a smaller fraction of the total emission budget, but a greater proportion of dissolved CH₄ is oxidized and thus the emitted CH₄ is more isotopically enriched. Our proposed conceptual model allows for these differences, where the higher efficiency of oxidation in larger rivers is enough to counter the higher prevalence of ebullition (Sawakuchi et al., 2021). Further, hyporheic flow is relatively less important in large rivers (Battin et al., 2008), potentially resulting in a greater proportion of CH₄ accumulating in bubbles in river sediments compared to streams. Even fewer studies exist with which to compare δD-CH₄ values. The streams in this study (−260‰) are δD-CH₄ enriched relative to Swedish lentic ecosystems (−310‰; Wik et al., 2020), which further supports the preferred emissions of oxidized CH₄ from streams. More closely comparing lotic δ¹³C-CH₄ signatures to wetlands of similar latitudes further supports relative isotopic enrichment in lotic systems. Latitudinal patterns have been observed in aquatic δ¹³C-CH₄ signatures, with enriched values typically characteristic of more tropical environments (Ganesan et al., 2018). The values measured in the present study are enriched relative to wetlands at a similar latitude (approx. −67‰; Fisher et al., 2017). Similarly, δ¹³C-CH₄ values measured in fluvial portions of the Amazon River basin (Sawakuchi et al., 2016; Sawakuchi et al., 2021) are enriched relative to other tropical wetlands (approx. −60‰; Brownlow et al., 2017). There remains uncertainty in how much aquatic CH₄ emissions are contributing to variation in atmospheric δ¹³C-CH₄ values (Turner et al., 2019), thus expanding measurement of aquatic CH₄ isotopes remains a valuable endeavor. While many additional measurements across ecosystems will be needed to confirm our proposed model, the possibility of a unique isotopic signature for streams and rivers is significant.