Water column respiration in the Yakima River basin is explained by temperature, nutrients, and suspended solids

Understanding aquatic ecosystem metabolism involves the study of two key processes: carbon fixation via primary production and organic C mineralization as total ecosystem respiration (ER_tot). In streams and rivers, ER_tot includes respiration in the water column (ER_wc) and in the sediments (ER_sed). While literature surveys suggest that ER_sed is often a dominant contributor to ER_tot, recent studies indicate that the relative influence of sediment-associated processes versus water column processes can fluctuate along the river continuum. Still, a comprehensive understanding of the factors contributing to these shifts within basins and across stream orders is needed. Here, we contribute to this need by measuring ER_wc and aqueous chemistry across 47 sites in the Yakima River basin, Washington, USA. We find that ER_wc rates vary throughout the basin during baseflow conditions, ranging from 0 to −7.38 g O₂ m⁻³ d⁻¹, and encompass the entire range of ER_wc rates from previous work. Additionally, by comparing to ER_tot estimates for rivers across the contiguous United States, we suggest that the contribution of ER_wc rates to reach-scale ER_tot rates across the Yakima River basin is likely highly variable, but we do not test this directly. We observe that ER_wc is locally controlled by temperature, dissolved organic carbon, total dissolved nitrogen, and total suspended solids, which explains 49 % of ER_wc variability across the basin using Least Absolute Shrinkage and Selection Operator (LASSO) regression. Our findings highlight the potential relevance of water column processes in aquatic ecosystem metabolism across the entire stream network and that these influences are likely not predictable simply by knowing the position in the stream network. Our results are generally congruent with previous work in terms of locally influential variables, suggesting that the observed variability and suite of associated environmental factors influencing ER_wc are potentially transferable across basins.

1 Introduction

Metabolism in streams and rivers includes both gross primary production (GPP) and ecosystem respiration (ER_tot) as fundamental processes that shape energy dynamics and nutrient cycling in riverine systems (Bernhardt et al., 2018). GPP and ER_tot impact biogeochemical cycling through the fixation and subsequent breakdown and processing of carbon (C) in aquatic ecosystems (Allan et al., 2021; Genzoli and Hall, 2016; Hall, 2016; Hall and Hotchkiss, 2017; Reisinger et al., 2016). Riverine metabolism is modulated by various environmental features, including physical and biogeochemical factors. Physical parameters include discharge, flow regimes, flow extremes, light availability, and temperature (Bernhardt et al., 2022; Hensley et al., 2019; Jankowski and Schindler, 2019; Nakano et al., 2022). Biogeochemical influences include the availability, amount, and composition of C and other nutrients (Bertuzzo et al., 2022; Garayburu-Caruso et al., 2020a; Mulholland et al., 2008; Reisinger et al., 2021). Additionally, watershed characteristics such as stream size or drainage area, hydrologic connectivity, watershed geomorphology, land use, and land cover further affect these metabolic processes (Bernot et al., 2010; Demars, 2019; Finlay, 2011; Jankowski and Schindler, 2019).

Reach-scale ecosystem metabolism encompasses biogeochemical processes that occur in both the water column and benthic and hyporheic sediments (Hall and Hotchkiss, 2017). Historically, metabolism studies focused on headwater streams, which are characterized by relatively large contact areas between surface water and the benthic sediments (Alexander et al., 2007; Battin et al., 2008; Findlay, 1995; Gomez-Velez et al., 2015; Mulholland et al., 2008; Peterson et al., 2001). Recent advances in computing power and the increased availability of high-resolution sensor data (e.g., dissolved oxygen, temperature, and river depth) have expanded the scope of metabolism studies beyond single small streams, enabling researchers to investigate the relative contributions of ER_sed and water column respiration (ER_wc) to ER_tot across diverse stream networks and orders. These efforts show that the proportion of ER_tot derived from ER_sed varies greatly across different sites, contributing from 3 % to 96 % of ER_tot (Battin et al., 2003; Fuss and Smock, 1996; Gagne-Maynard et al., 2017; Jones, 1995; Kaplan and Newbold, 2000; Naegeli and Uehlinger, 1997). This observed variability in the fraction of ER_tot derived from ER_sed indicates that ER_wc may be important in certain places and at certain times.

Water column processes, including nutrient cycling, occur at considerable rates and become increasingly important as rivers grow in size, marking a transition from benthic-dominated to water-column-dominated processing (Gardner and Doyle, 2018; del Giorgio and Williams, 2005; Reisinger et al., 2015, 2016), and they may influence ecosystem respiration such that we would expect faster ER_wc with greater GPP due to increases in C (Hall et al., 2016; Mejia et al., 2019). Additionally, greater N processing in the water column with increasing stream order (Wang et al., 2022) may suggest that water column biogeochemical processing increases along the stream network. Despite these trends, even as rivers increase in size, the relative contribution of ER_wc to ER_tot remains variable, likely in response to changing environmental conditions (Genzoli and Hall, 2016; Reisinger et al., 2021; Ward et al., 2018). This highlights a key knowledge gap: while the role of the water column in reach-scale processes such as GPP and ER_tot likely fluctuates along the river network, this relationship remains poorly understood.

We contribute to addressing this knowledge gap by investigating the spatial variation of ER_wc in the Yakima River basin, Washington, USA. The Yakima River basin is representative of the Columbia River basin, one of the largest river basins in the contiguous United States (CONUS), which spans the northwest region of CONUS. The Yakima River basin encompasses climatic regimes, biomes, physical settings, and land use conditions commonly found throughout the Columbia River basin and the western CONUS. Using the environmental diversity of the Yakima River basin, our goal is to generate knowledge of ER_wc that could be transferable across the Columbia River basin and potentially beyond. We focus on ER_wc during summer baseflow conditions and specifically (1) compare ER_wc from the Yakima River basin to published ER_wc and ER_tot from other systems; (2) test the hypothesis that ER_wc will be faster as one moves down the stream network; and (3) compare variables that explain the variation in ER_wc to those found as explanatory in previous studies. To address these objectives, we estimate ER_wc and measured surface water chemistry at 47 sites across the Yakima River basin during the summer of 2021. Our estimates of ER_wc span all previously reported rates, and the most important explanatory variables align with previous studies in a way that indicates predominance of local controls.

2 Methods

2.1 Methods overview

Field sites in the Yakima River basin were selected to be representative of biophysical attributes of the larger Columbia River basin. For this, we grouped all catchments in the Columbia River basin into six classes sharing similar landscape characteristics using key biophysical attributes and selected sites in the Yakima River basin from each of the six classes. The final field locations spanned six Strahler stream orders and a wide range of land cover types and physical settings. We used dark bottle incubations and collected surface water chemistry samples to study the spatial variability of ER_wc at the basin scale with respect to environmental conditions during summer baseflow conditions in 2021. We also compared ER_wc observed in the Yakima River basin against literature ER_wc and ER_tot values to understand how the Yakima River basin relates to streams and rivers across the world. We used Least Absolute Shrinkage and Selection Operator (LASSO) regression to evaluate the relationship between ER_wc and drainage area, stream temperature, surface water chemistry, and organic matter putative biochemical transformations as a proxy for the diversity of reactions occurring in upstream reaches to determine the primary factors influencing ER_wc throughout the Yakima River basin. All analyses were performed using R Statistical Software (v4.2.0) (R Core Team, 2022). All data generated from the sampling study, including data not evaluated in this paper, are publicly available.

2.2 Watershed characterization and site selection

The Yakima River basin is the fifth-largest basin in the Columbia River basin and is located entirely within the state of Washington, USA. The basin is roughly 16 000 km² and spans forested mountainous regions in the west to arid valleys and plains in the east. The basin has a diversity of land covers and land uses dominated by shrubland, forest, and agriculture. Annual precipitation ranges from up to 350 cm in the west to 25 cm in the east (Vano et al., 2010).

To enable further testing of the transferability of the study results to catchments throughout the Columbia River basin, we strategically selected sampling sites in the Yakima River basin based on their biophysical (e.g., hydrology, topography, vegetation type) characteristics. This was done by first grouping all National Hydrography Dataset Plus Version 2.1 (NHDPlusV2.1) catchments (U.S. Geological Survey, 2019a) in the Columbia River basin (n=181 531) into six classes sharing similar landscape characteristics using cluster analysis. To capture the variability in biophysical settings found across the Columbia River basin, we selected 16 key attributes as input variables to the cluster analysis, including climate, vegetation structure and function, topography, and wildfire potential (Table S1 in the Supplement). We then selected multiple sites within each of the six Columbia River basin classes. Existing, readily available geospatial data came from multiple sources, including the NASA Moderate Resolution Imaging Spectroradiometer (eMODIS) Remote Sensing Phenological (RSP) data (U.S. Geological Survey, 2019b); NASA MODIS land cover type (Friedl and Sulla-Menashe, 2019); NASA MODIS normalized difference vegetation index (NDVI), fraction of photosynthetically active radiation (FPAR, %), and leaf area index (LAI, m² m⁻²) (Myneni et al., 2015); NASA MODIS total evapotranspiration (ET, kg H₂O m² d¹) (Running et al., 2017); NASA MODIS terrestrial net primary productivity (NPP, kg C m⁻² yr⁻¹) and terrestrial net ecosystem productivity (NEP, kg C m⁻² yr⁻¹) data (Running and Zhao, 2019); PRISM precipitation data (PRISM Climate Group, Oregon State University, 2014); NHDPlusV2.1 stream length and catchment boundaries (U.S. Geological Survey, 2019a); USGS National Elevation Dataset (NED) $\mathrm{1}/\mathrm{3}$ arcsec Digital Elevation Model topography data (National Elevation Dataset (NED), 2023); USFS Wildfire Hazard Potential (WHP) data (Dillon, 2018); and Landscape Fire and Resource Management Planning Tools (LANDFIRE) existing vegetation percent cover (%) and height (m) data (Dillon and Gilbertson-Day, 2020).

We used a k-means clustering algorithm using the kmeans function within the “stats” package in base R to group NHDPlusV2.1 catchments with similar properties using the normalized, statistical moments (minimum, maximum, mean, and standard deviation (SD)) for 70 geospatial variables within each NHDPlusV2.1 catchment (Table S1 in the Supplement) as input. To calculate statistical moments for each variable, we summarized geospatial data types at the NHDPlusV2.1 catchment level using two different methods: zonal statistics for continuous raster data and tabulation for vector data. Zonal statistics calculate statistical moments by individual catchment polygon. Tabulation calculates total length or area of a particular vector feature within each individual catchment polygon. We evaluated 13 different sets of variable–statistical moment combinations for use in the cluster analysis and selected variable set 8, which included the zonal mean and zonal standard deviation for 70 variables (n=140) (Table S2). Once the data for variable set 8 were summarized at the NHDPlusV2.1 catchment level, we calculated z-scores (z) for each geospatial variable. Resultant z-scores for variable set 8 were fed into the k-means classifier, which iteratively adds each catchment to one of n clusters, with n being set by the user (n=15 in this study), using the Euclidean distance to minimize the within-cluster distance and maximize the between-cluster distance. We ran multiple iterations of the cluster analysis using 2–15 clusters with the mean and standard deviation of all variables. To visualize the reduction in within-cluster variation between iterations 1–15, we generated elbow plots by plotting the within-cluster sum of squares (WCSS) value against the total number of catchments in a cluster and selected six clusters as the suitable number of clusters that minimized map visual complexity enough to guide manual site selection while maintaining a level of variation in key biophysical characteristics representative of the Columbia River basin. Clusters 1 and 3–6 were categorized according to tree height, precipitation, and elevation (Tables 1 and S3). Cluster 2 was categorized as “Water dominated” and was not used for selecting sites. The cluster analysis results were then used to guide the selection of 47 field sites distributed across Strahler stream orders 2–7 (the highest-order stream in the Yakima River basin) that spanned the basin and captured the variation in biophysical characteristics represented by clusters 1 and 3–6 (Fig. S1 in the Supplement). First-order and other non-perennial streams were not sampled due to the lack of flow during summer baseflow or because the baseflows were too low to support sampling. We attempted to include logistical considerations in the model-based site selection, but this task proved impractical, and field-scouting trips were needed to refine site selections. Day-of-sampling changes to the sampling plan were made on-the-fly when the Schneider Springs Fire started at the Okanogan–Wenatchee National Forest. Fire activity and road closures restricted access to a large portion of the Yakima River basin, primarily in the Tieton River and American River watersheds located in the midwestern portion of the basin.

Table 1 Cluster analysis results characterizing NHDPlusV2.1 catchments across the Columbia River basin and Yakima River basin with similar biophysical and hydrologic characteristics, and the number and percentage of sites in each basin.

“CRB drainage area” is the percentage of the total drainage area of the Columbia River basin that was classified in each cluster. “YRB drainage area” is the percentage of the total drainage area of the Yakima River basin that was classified in each cluster. “YRB sites per cluster” is the total number of field sites in the Yakima River basin (n=47) located in each cluster. “Percent YRB sites per cluster” is the percentage of the total number of sampling sites in the Yakima River basin located in each cluster.

2.3 Water column respiration data collection

We measured ER_wc (g O₂ m⁻³ d⁻¹) in triplicate for 2 h at each site between 30 August and 15 September 2021 using a modified “semi-in situ” dark bottle incubation (Genzoli and Hall, 2016) (Fig. 1a). Calibrated dissolved oxygen (DO) sensors (miniDOT Logger; Precision Measurement Engineering, Inc.; Vista, CA, USA) recorded DO concentration (mg L⁻¹) and temperature (°C) at 1 min intervals in 2 L dark bottles (Nalgene™ Rectangular Amber HDPE bottles; Thermo Fisher Scientific; Waltham, Massachusetts, USA) (Fulton et al., 2022). Bottle necks were slightly widened (1 to 2 mm) to accommodate the diameter of the DO sensor.

At the start of each sampling day, the DO sensors and all sampling equipment were placed in a cooler with blue ice packs to keep them cool and minimize the time needed at each site for the sensors to equilibrate with the similarly cool river water temperatures. Upon arrival at each site, the bottles were rinsed three times with river water and then filled by wading as close to the thalweg as possible, submerging the bottles below the river surface, and rolling them 360° while held upright underwater to ensure no air bubbles were present in the bottles (Fig. 1a). The bottles were secured upright in a cooler filled with river water, placed in the shade on the streambank, and allowed to equilibrate for 20 min. Following the 20 min equilibration period, the bottles were emptied and refilled with fresh river water and a small, battery-powered mixing device (Underwater Motor, Item Number 7350; Playmobil; Shanghai, China; rechargeable AA NiMH battery; Amazon; Seattle, Washington, USA), and the DO sensor was gently inserted (sensor face-up) in the bottles to minimize the trapping of air bubbles in the bottles. The bottles were capped underwater and returned to the water-filled cooler. The bottles were incubated for 2 h, and the river water surrounding the bottles in the cooler was replenished every 20 min to maintain in situ temperature.

Figure 1 Modified semi-in situ dark bottle incubation method and example study sites. (a) Underwater photograph of DO sensor being inserted into an incubation bottle filled with river water and the mixing device. Panels to the right emphasize the diversity of environmental settings covered in this study. (b) North Fork Teanaway River (site S19E), Kittitas County, Washington, September 2021. Site S19E is classified as a mesic, high-elevation site dominated by tree canopy (cluster 1; see Tables 1 and S3 and Fig. S1). (c) Yakima River at Mabton (site T02), Yakima County, Washington, September 2021. Site T02 is classified as a mesic, middle-elevation site dominated by tree canopy (cluster 5; see Tables 1 and S3 and Fig. S1).

2.4 Surface water chemistry sample collection and analysis

Filtered surface water samples were collected at each site for dissolved inorganic C (DIC, mg L⁻¹); dissolved organic C (DOC, mg L⁻¹); total dissolved N (TDN, mg L⁻¹); anions, including nitrate (NO${}_{\mathrm{3}}^{-}$, mg L⁻¹), chloride (Cl⁻, mg L⁻¹), and sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, mg L⁻¹); and DOM chemistry using a 50 mL syringe and 0.22 µm sterivex filter (MilliporeSigma™ Sterivex™ Sterile Pressure-Driven Devices; MilliporeSigma™; Burlington, Massachusetts, USA) (Grieger et al., 2022). Samples were collected in triplicate from 50 % of the water column depth. Prior to sample collection, filter assemblies were rinsed once by pushing 5 mL of river water through the filter. DIC, DOC, and TDN samples were filtered into 40 mL amber glass vials (Amber Clean Snap Vials; Thermo Fisher Scientific; Waltham, Massachusetts, USA). DIC samples were collected by attaching a sterile 18 g needle (BD General Use and PrecisionGlide Hypodermic Needles; Becton, Dickinson, and Company; Franklin Lakes, NJ, USA) to the sterivex filter and pushing three vial-volumes of river water (∼ 150 mL) slowly through the syringe to prevent the introduction of air bubbles to the sample, allowing the vials to overflow continuously. When the final 50 mL of river water was pushed through the syringe, the vials were capped with a surface tension dome of water to ensure no headspace. Samples collected for ion analysis were filtered into a 15 mL conical tube (Olympus™ Plastics; Genesee Scientific; Morrisville, NC, USA). Samples collected for DOM chemistry were filtered into pre-acidified (85 % phosphoric acid, H₃PO₄) 40 mL amber vials (Amber Clean Snap Vials; Thermo Fisher Scientific; Waltham, Massachusetts, USA) (Grieger et al., 2022). One unfiltered grab sample for total suspended solids (TSS, mg L⁻¹) was collected using a pre-washed 2 L amber bottle (Nalgene™ Rectangular Amber HDPE Bottles; Thermo Fisher Scientific; Waltham, Massachusetts, USA). TSS bottles were rinsed three times with river water prior to sample collection. All samples were stored on ice in the field and then refrigerated at 4 °C before shipping for analysis to the Pacific Northwest National Laboratory (PNNL) Marine and Coastal Research Laboratory in Sequim, Washington (DOC and DIC), and PNNL Biological Sciences Facility Laboratory in Richland, Washington (TSS, ions, and DOM). The TSS samples were analyzed within 1 week of collection, DOC and TDN were measured within 2 weeks of collection, DIC was measured within 1 month of collection, and ion and DOM samples were frozen (−20 °C) upon receiving until analysis.

DOC, TDN, and DIC were measured on a Shimadzu TOC-L total organic carbon analyzer. DOC was measured as non-purgeable organic C (NPOC). Anion concentrations were determined quantitatively on a Dionex ICS-2000 anion chromatograph with an AS40 autosampler using one replicate. An isocratic method was used with 23 mM KOH eluent at 1 mL min⁻¹ at 30 °C. The analytical column was an IonPac AS18 (4×250 mm, Dionex catalog # 060549). The suppressor was an ADRS 600 set at 57 mA (4 mm, self-regenerating, Dionex catalog # 088666). Concentrations below the limit of detection (LOD) of the instrument, or below the standard curve, were flagged (Grieger et al., 2022). For other samples below the lowest standard value (TDN: 0.1 mg L⁻¹, NO${}_{\mathrm{3}}^{-}$: 0.07 mg L⁻¹), one-half of the lowest standard value was used (TDN: 0.05 mg L⁻¹, NO${}_{\mathrm{3}}^{-}$: 0.035 mg L⁻¹) for statistical analysis. For samples below the limit of detection (TDN LOD: 0.07 mg L⁻¹; NO${}_{\mathrm{3}}^{-}$ LOD: 0.07 mg L⁻¹) but above the lowest standard, one-half of the LOD value (TDN: 0.035 mg L⁻¹; NO${}_{\mathrm{3}}^{-}$: 0.035 mg L⁻¹) was used for analysis. Phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$) was measured; however, over two-thirds of samples showed values below detection, and thus the analyte was not used in subsequent analyses. Pairwise differences between the NPOC, TDN, and DIC measurements from all replicates were calculated. The sample that had the largest difference from the other samples was removed if the coefficient of variation was greater than 30 %. This coefficient of variation threshold for sample removal is based on inspecting histograms of these data types and determining the point at which sites likely contain anomalous outlier values. Parameter mean values for each site were then calculated from the remaining replicates.

TSS samples were filtered in the laboratory through a pre-weighed and pre-combusted 4.7 cm, 0.7 µm GF/F glass microfiber filter (Whatman™ glass microfiber filters, Grade 934-AH^®; MilliporeSigma; Burlington, Massachusetts, USA). After water filtration, the filter and filtration apparatus were rinsed with 30 mL of ultrapure Milli-Q water (Milli-Q^® IQ Water Purification System; MilliporeSigma; Burlington, Massachusetts, USA) to ensure that all residue was captured by the filter. The filter was placed in foil and oven-dried overnight at 45 °C. TSS (mg L⁻¹) was calculated as the difference between the weight (mg) of the filter before and after filtration of the water sample divided by the volume of water filtered (L). For samples below the LOD, one-half of the LOD value (LOD: 0.24 mg L⁻¹) was used for analysis.

2.5 DOM chemistry via ultra-high-resolution mass spectrometry and biochemical transformations

Organic matter chemistry was characterized via ultra-high-resolution mass spectrometry using a 12 Tesla (12T) Bruker solariX Fourier transform ion cyclotron resonance mass spectrometer (FTICR-MS) at the PNNL Environmental Molecular Sciences Laboratory in Richland, Washington, following methods described in Garayburu-Caruso et al. (2020b). Measured DOC concentrations were used to normalize the DOC concentration of the sample to 1.5 mg C L⁻¹ prior to further processing. Samples were thawed in the dark at 4 °C overnight before acidifying to pH 2 using 85 % H₃PO₄. Samples were then subjected to solid phase extraction (SPE) using Bond Elut PPL cartridges (Agilent; Santa Clara, CA, USA) following protocols employed by Dittmar et al. (2008). Extracted samples were run in the FTICR-MS with a standard electrospray ionization source in negative mode. Data were collected with an ion accumulation time of 0.08 s. Bruker Daltonik Data Analysis version 4.2 was used to convert raw spectra to a list of molecular compound mass-to-charge ratios ($m/z$) with a signal-to-noise ratio ($S/N$) threshold set to 7 and absolute intensity threshold to the default value of 100. Peaks were aligned (0.5 ppm threshold), and molecular formulas were assigned using the Formultitude (previously named Formularity) software, with $S/N>\mathrm{7}$ and mass measurement error < 0.5 ppm (Tolić et al., 2017). The compound identification algorithm takes into consideration the presence of C, H, O, N, S, and P and excludes other elements. The aligned and calibrated data were further processed using ftmsRanalysis (Bramer et al., 2020). Replicate samples were merged into one site, where peaks in a sample were retained if they were present in at least one of the replicates. DOM biochemical transformations were inferred following methods previously employed by Danczak et al. (2023), Fudyma et al. (2021), Garayburu-Caruso et al. (2020a), Ryan et al. (2024), and Stegen et al. (2018). In summary, we calculated pairwise mass differences between every peak in a sample regardless of the molecular formula assigned and compared that mass difference to a list of 1255 molecular masses associated with commonly observed biochemical transformations (Table S4). Biochemical transformations allow you to infer the number of times the mass that corresponds to a specific molecule is gained or lost. For example, if a mass difference between two peaks corresponded to 128.095, that would correlate to the loss or gain of the amino acid lysine (see Table S4). We further calculated the total number of DOM transformations per site and the total number of DOM transformations normalized by the number of peaks present in the site (i.e., “normalized DOM transformations”).

2.6 DO sensor data cleaning, processing, and analysis

We extracted the raw DO concentration (mg O₂ L⁻¹) and temperature (°C) sensor data for each site and plotted DO and temperature against incubation time for each set of triplicate incubations (n=141). The plots were visually inspected to (a) confirm that the temperature sensors were at equilibrium with the river temperature when the 2 h incubation test period began and (b) identify data gaps, outliers, and other data anomalies. Following the visual inspection of plots, the first 5 min of the time series was removed, and then the data were trimmed to 90 min to account for anomalies due to emptying and refreshing river water in the bottles and to ensure that all sites had the same incubation time. The sensor data distributions were also evaluated using violin plots for each site.

ER_wc rates for individual triplicate incubation samples were calculated as the slope of the linear regression between the DO sensor data and the incubation time, which was converted to units of g O₂ m⁻³ d⁻¹. All samples met the normalized root mean square error (NRMSE) criteria of ≤ 0.01 (Shcherbakov et al., 2013). Mean ER_wc for each site and the global mean and variance were then calculated from the samples (n=141). Nearly one-fifth of the ER_wc values were slightly positive. Positive respiration rates are biologically unrealistic; however, positive values less than 0.5 g O₂ m⁻³ d⁻¹ are difficult to distinguish from 0 (Appling et al., 2018b). Thus, we changed positive ER_wc values less than 0.5 g O₂ m⁻³ d⁻¹ to 0 for analysis and removed values greater than 0.5 g O₂ m⁻³ d⁻¹ (n=2). ER_wc values greater than 0.5 g O₂ m⁻³ d⁻¹ were observed when the DO concentration in the bottle started near 5 mg O₂ L⁻¹ and increased over the 2 h incubation period. The increase in concentration and the high, positive respiration rate are likely due to the diffusion of DO through the bottle walls. The ER_wc values are reported in volumetric units (g m⁻³ d⁻¹) as opposed to areal units (g m⁻² d⁻¹) due to difficulties in obtaining high-quality depth data across all field sites, spanning small headwater streams to large mainstem rivers.

2.7 Relationship of water column respiration rates to watershed characteristics and surface water chemistry

We evaluated the relationship between ER_wc, watershed characteristics, physical parameters, and surface water chemistry using LASSO regression models, which perform variable selection and model regularization, to establish the suite of explanatory variables that most influence variation in ER_wc across the Yakima River basin. We observed that several model input variables had skewed distributions; thus, a cube root transformation was applied to all variables to reduce the impact of high leverage points in the regression analysis. Further, all data were standardized as z-scores before analysis to ensure that all the data were in the same quantitative range. The β coefficients reported for each variable were calculated by performing LASSO regression using the glmnet function in R (Friedman et al., 2010) over 100 random seeds, normalizing to the maximum β coefficient in each regression, and averaging the normalized β coefficients across the 100 iterations. The minimum penalty parameter (λ) determined by cross validation was used in each regression. Because LASSO regression was used for exploratory analysis, not prediction, the data were not split into training and testing sets. LASSO does not inherently estimate R², so we calculated it using the total sum of squares and residual sum of squares for each fitted model, as traditionally done with standard multiple regression. For estimation of the residual sum of squares, we used predicted values of ER_wc based on the explanatory variables used in the model. The R² estimates were used to estimate how much variation in ER_wc was explained by each of the LASSO models. The standard deviations of the β coefficients were compared to the mean values of the β coefficients to confirm that the most important variables were relatively consistent across seeds. Total drainage area (km²) was defined as the total upstream drainage area from each site and was extracted for each site from the NHDPlusV2.1 stream database using site latitude and longitude. The stream order for each site was extracted as the reach attribute “StreamOrde” from the NHDPlusV2.1 stream database, which is a modified version of Strahler stream order (Blodgett and Johnson, 2022; U.S. Geological Survey, 2019a; Willi and Ross, 2023). To evaluate whether the directionality of relationships observed in the LASSO regression was consistent with univariate relationships, we used Pearson correlations between ER_wc, drainage area, water chemistry, and environmental factors; these correlations were calculated using the cor function in R.

2.8 Comparison to published water column respiration rates

To contextualize the magnitude of observed ER_wc rates in the Yakima River basin, we compared our results to published literature values of ER_wc (n=118) (Table S5) and ER_tot (n=208). Published ER_wc values were converted to volumetric units (g O₂ m⁻³ d⁻¹) using standard unit conversions. For example, molar values (µmol O₂ L⁻¹ h⁻¹) as in Devol et al. (1995) and Quay et al. (1995) were corrected using the molar mass of oxygen and standard time conversions. When ER_wc was reported with respect to C or carbon dioxide (CO₂), as in Ellis et al. (2012) and Ward et al. (2018), the conversions provided in the text were used to convert to an O₂ basis. Areal estimates of ER_wc (g O₂ m⁻² d⁻¹), as in Genzoli and Hall (2016) and Reisinger et al. (2021), were converted to volumetric units by multiplying by 1/depth (m⁻¹) using same-day depth data for each reach studied. We also compared our ER_wc values to daily reach-averaged estimates of ER_tot (n=490 907) for 356 rivers and streams across the CONUS by using the datasets published in Appling et al. (2018a, b) and Bernhardt et al. (2022), where ER_tot was estimated by a single-station, open-channel approach using the streamMetabolizer package in R (Appling et al., 2018a, b). For our comparative analysis, we used the cleaned, gap-filled data from Bernhardt et al. (2022) (n=208). The Bernhardt et al. (2022) sites are a subset of the Appling et al. (2018a, b) dataset generated through a robust data quality analysis to remove sites potentially affected by process or observation error. For comparison with our ER_wc values, we first averaged the Bernhardt et al. (2022) ER_tot areal units (g O₂ m⁻² d⁻¹) at each site. Then, the average ER_tot values were converted to volumetric units by calculating the average river depth per site from the Appling et al. (2018a, b) dataset and multiplying the average ER_tot by 1/depth.

3 Results and discussion

3.1 Yakima River basin ER_wc rates spanned literature values

At baseflow conditions, ER_wc varied widely across the Yakima River basin. The linear regression models for each triplicate set of DO sensor measurements were well-fit to the data, and all sites met the criteria for NRMSE ≤ 0.01 (Figs. S2, S3). We observed consistency across triplicate measurements, illustrating that the method was effective in providing repeatable estimates of ER_wc rates throughout the Yakima River basin (Figs. S2, S3). After removing positive respiration rates > 0.5 g O₂ m⁻³ d⁻¹, which were associated with diffusion effects on DO, and turning small positive rates to 0, the ER_wc rates ranged from −7.38 to 0 g O₂ m⁻³ d⁻¹, with a median value of −0.58 g O₂ m⁻³ d⁻¹ (mean: −0.84 g O₂ m⁻³ d⁻¹, standard deviation = 1.23 g O₂ m⁻³ d⁻¹) (Fig. 2a).

The values of ER_wc observed in our study spanned the range of published literature values (Fig. 2; Table S5). From 118 published measurements of ER_wc across the CONUS and the Amazon River basin, ER_wc ranged from −4.63 to −0.07 g O₂ m⁻³ d⁻¹. We compared median values, rather than means, across studies because medians are more appropriate for skewed distributions and are less sensitive to outliers in the data. The median ER_wc from this study (−0.58 g O₂ m⁻³ d⁻¹) is slower than the median of literature-reported ER_wc values (−0.96 g O₂ m⁻³ d⁻¹). However, the fastest ER_wc rate in the Yakima River basin (−7.38 g O₂ m⁻³ d⁻¹) exceeded the fastest reported literature value (−4.63 g O₂ m⁻³ d⁻¹) (Reisinger et al., 2021). Reisinger et al. (2021) measured ER_wc in 15 mid-sized rivers across basins with differing turbidity levels and nutrient concentrations, finding a similar median ER_wc (−0.60 g O₂ m⁻³ d⁻¹) to that in this study. In the Klamath River, the median ER_wc (−0.51 g O₂ m⁻³ d⁻¹) was also similar to that in the Yakima River basin. However, ER_wc doubled following summer cyanobacteria blooms, emphasizing the temporal variability in water column processes with changing environmental conditions (Genzoli and Hall, 2016). In the Amazon basin, literature comparisons varied, with median ER_wc measurements similar to those found in the Yakima River basin in some studies (Devol et al., 1995; Ellis et al., 2012; Quay et al., 1995) and faster than those in the Yakima River basin in others (Ward et al., 2018). Ward et al. (2018) highlighted the importance of mixing in large rivers, noting that previous measurements of aquatic respiration in large tropical rivers, such as those measured in Quay et al. (1995) and Devol et al. (1995), may underestimate microbial respiration contribution due to lack of mixing during rate measurements. While comparisons across study medians are variable, the observation that ER_wc in the Yakima River basin spans – and exceeds – reported literature values highlights the potential for using it as a test basin for understanding and uncovering transferable principles linked to stream metabolism.

While ER_tot estimates are not available across the Yakima River basin at the time of ER_wc estimation for this paper, the measured ER_wc rates spanned a large fraction of CONUS-scale ER_tot rates estimated by Appling et al. (2018a, b) and Bernhardt et al. (2022). ER_tot rates are reach-scale estimates of stream metabolism derived from time series measurements of DO. This method assumes well-mixed conditions such that sensor measurements represent homogenous reach observations. Under well-mixed conditions, ER_wc measurements from dark bottle incubations are also representative of reach-scale processes (Genzoli and Hall, 2016). The median ER_tot for 208 CONUS measurements was −5.25 g O₂ m⁻³ d⁻¹ with a range from −36.55 to −3.73 g O₂ m⁻³ d⁻¹. The median ER_wc rate (−0.58 g O₂ m⁻³ d⁻¹) observed in the Yakima River basin was 11 % of the median ER_tot (Fig. 2). The fastest ER_wc rate in the Yakima River basin (−7.38 g O₂ m⁻³ d⁻¹) was faster than the median ER_tot (Fig. 2). While both ER_tot and ER_wc measurements span a range of stream conditions, we acknowledge that we did not compare these rates directly at the same places and times. However, given the overlap of ER_wc from the Yakima River basin with CONUS-scale ER_tot, we suggest that ER_wc could typically represent a small fraction of ER_tot but may occasionally have larger contributions across the Yakima River basin. If we had observed consistently very slow ER_wc across the Yakima River basin, there would be little overlap with the literature ER_tot values, and we would have inferred consistently small contributions of ER_wc to ER_tot. In comparison, Genzoli and Hall (2016) observed that before summer cyanobacteria blooms, ER_wc contributed around 10 % of ER_tot in sites along the Klamath River, with the contribution of ER_wc to ER_tot increasing following cyanobacteria blooms. Additionally, Reisinger et al. (2021) found that ER_wc was not the dominant contributor to ER_tot in mid-sized rivers, except at sites with low ER_tot (mean ER_wc contributions to ER_tot: 35 %, range 2 %–81 %). While these studies have shown spatiotemporal variability of the contributions of ER_wc to ER_tot, exploring these relationships in the Yakima River basin requires further research where ER_tot is measured in conjunction with ER_wc.

Figure 2 Water column respiration data from the Yakima River basin (ER_wc (this study); n=45), published water column respiration rates (ER_wc (Lit); n=118), and reach-scale estimates of ecosystem respiration by Appling et al. (2018a, b) and Bernhardt et al. (2022) (ER_tot; n=208). (a) Kernel density plots of ER_wc from the Yakima River basin (this study), published ER_wc rates (Lit) that have been converted to the same units as in this study (g O₂ m⁻³ d⁻¹), and published reach-scale ER_tot (Lit) from Bernhardt et al. (2022) that have been converted to volumetric units using depth data from Appling et al. (2018a). The left y axis is for ER_wc values. The right y axis is for ER_tot values. The vertical blue line is the median ER_wc observed in the Yakima River basin (−0.58 g O₂ m⁻³ d⁻¹). The vertical red line is the median ER_wc values from studies in rivers across the CONUS and the Amazon River basin (−0.96 g O₂ m⁻³ d⁻¹). The vertical black line is the median ER_tot value (−5.25 g O₂ m⁻³ d⁻¹). (b) Boxplots of published ER_wc and ER_wc from the Yakima River basin. The blue horizontal dashed line represents median ER_wc in the Yakima River basin. The red horizontal dashed line represents median ER_wc from literature values.

3.2 Water column respiration rates varied weakly with drainage area and stream order

We observed a relatively weak correlation between ER_wc and drainage area across the Yakima River basin, though it does suggest that ER_wc becomes faster while moving down the stream network (Fig. 3). In latter sections, we use multivariate analysis for further evaluation of the relationships between ER_wc and explanatory variables, effectively competing drainage area against other variables as a way to more deeply evaluate our hypothesis. The lack of a strong connection between ER_wc and drainage area is somewhat surprising, as a strong relationship could emerge from downstream C transport as well as increasing autochthonous C inputs due to increasing temperature and light availability, providing an additional substrate for microbial respiration (Finlay, 2011; Webster, 2007). The fastest observed ER_wc rate in the Yakima River basin occurred in an agriculturally influenced, low gradient, fifth-order stream, as opposed to our hypothesis of ER_wc being fastest in the highest stream orders (Fig. 3). The conditions at this sampling location were not representative of the whole drainage area, as areas upstream of this site are mountainous with little human influence. This finding suggests that localized factors, not upstream conditions or drainage area, provide primary controls over ER_wc. Anthropogenic impacts, such as from agriculture and urbanization, can alter nutrient dynamics and flow regimes in these areas, influencing biogeochemical processes such as ER_wc (Bernot et al., 2010). Additionally, while we report ER_wc on a volumetric basis, we acknowledge that this approach does not account for variation in water column depth along the river continuum. As river depth increases downstream, we expect the areal contribution of water column processes to also increase because areal contributions integrate across the whole water column (Wang et al., 2022). The relatively weak correlation between volumetric-based ER_wc and drainage area in the Yakima River basin likely reflects the interplay of multiple factors, including spatially variable local conditions, underscoring the complex controls on ecosystem processes in this region.

Figure 3 ER_wc across the Yakima River basin and its relationship with total drainage area. (a) Map of land use/land cover classes in the Yakima River basin with ER_wc values (g O₂ m⁻³ d⁻¹) overlaid. Faster rates are indicated by larger circle diameters. The fastest rate is indicated by the yellow circle. The map was generated using the Free and Open Source QGIS (v. 3.16.1 and v. 3.26.0). Map data include catchment boundaries and hydrography from the National Hydrography Dataset Plus (NHDPlusV2.1) (U.S. Geological Survey, 2019a) and 2016 land use/land cover data from the National Land Cover Dataset (U.S. Geological Survey, 2024). (b) Scatter plot of cube-root-transformed ER_wc relative to cube-root-transformed total drainage area, with points colored by stream order. The Pearson correlation coefficient (r) is provided on the panel.

3.3 Higher temperatures and nutrient concentrations are associated with faster ER_wc

Regression analyses showed that ER_wc in the Yakima River basin varied with chemical and physical water quality parameters. TDN, temperature, DOC, and TSS emerged as key variables in the LASSO regression, whereby ER_wc was faster with higher values of all these variables (Table 2). The LASSO regression explained 49 % of the variation in ER_wc (Table 2). LASSO results are similar to univariate relationships, whereby DOC, TDN, temperature, and TSS had the strongest correlations with ER_wc ($r=-\mathrm{0.46}$ to −0.63) (Figs. 4, S4) and all correlations were qualitatively in the same direction, as indicated by the LASSO β coefficients. Changing positive ER_wc values less than 0.5 g O₂ m⁻³ d⁻¹ to 0 did not change the overall interpretation of univariate or multivariate relationships (Figs. S4, S5, Table S4). Collectively, the relative importance of these variables suggests that ER_wc is not controlled by a single variable, and instead multiple factors (i.e., nutrient concentrations, suspended particles, and temperature) are simultaneously linked to ER_wc.

Collinearity between LASSO variables could result in one variable being retained in the LASSO model over another. We used LASSO regressions across 100 random seeds, averaging the model coefficients, to help minimize spurious outcomes. This revealed relatively small standard deviations of β coefficients compared to mean β coefficient values, indicating that the four most important variables are consistent across seeds, even when one variable is chosen over another (Table 2). For example, total drainage area was correlated with nutrient concentrations and temperature (Fig. S4), which were retained as more directly explaining variation in ER_wc in the LASSO regression. Additionally, while total drainage area showed a negative univariate correlation with ER_wc (Fig. 3b), it showed a slight positive relationship with ER_wc in the LASSO regression. This suggests that total drainage area likely acts as a proxy for regional watershed processes that influence ER_wc directly, like nutrients and temperature, rather than a causal relationship (Caissie, 2006; Manning et al., 2020). Similarly, TDN was strongly correlated with other explanatory variables, such as NO${}_{\mathrm{3}}^{-}$, Cl⁻, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$, likely reflecting an increase in agricultural inputs that, in turn, lead to faster ER_wc through supporting microbial metabolism (Bernot et al., 2010). Including phosphorus data could further improve model performance, as phosphorus is often a limiting factor for microbial growth in freshwater rivers (Carroll, 2022). Phosphorus limitation is likely in the Yakima River basin, as more than two-thirds of the phosphorus concentrations were below instrument detection, leading to its exclusion from the analysis. These results underscore the importance of interpreting LASSO results within the context of all explanatory variables used in the model, particularly in large, heterogenous catchments.

Table 2 β coefficients from LASSO analyses for explaining ER_wc across the Yakima River basin. ER_wc and all explanatory variables were cube-root-transformed and standardized as z-scores. LASSO was performed over 100 seeds, and β coefficients for each variable were normalized to the maximum β coefficient in each seed and averaged across all seeds for the reported values. Values of 0 indicate that, while the variable was included in the model, it was deemed not influential in predicting model outcomes and thus was not assigned a β coefficient.

Figure 4 Scatter plots of cube-root-transformed variables that were important in the LASSO regression. Cube-root-transformed ER_wc is the y axis for all panels. (a) Cube-root-transformed total dissolved nitrogen (TDN); (b) cube-root-transformed temperature; (c) cube-root-transformed dissolved organic carbon (DOC); and (d) cube-root-transformed total suspended solids (TSS). Pearson correlation coefficients (r) are provided on each panel.

Faster ER_wc with increasing TDN, temperature, DOC, and TSS in the Yakima River basin is expected, as nutrients and temperature can impact variation in stream metabolism (Ardón et al., 2021; Bernot et al., 2010; Honious et al., 2021; Hornbach, 2021; Nakano et al., 2022). In-stream metabolism relies on terrestrially derived and internally fixed inputs of DOC, which supports heterotrophic metabolism that degrades and removes organic C inputs through respiration (Hall et al., 2016; Hotchkiss and Hall, 2014; Plont et al., 2022). Faster ER_tot and ER_wc have been reported with increases in DOC (Bernot et al., 2010; Ellis et al., 2012). However, elevated DOC does not always correspond to greater ER_tot, as discharge and residence time also affect C dynamics (Ulseth et al., 2018). In addition to DOC, suspended sediment can regulate ecosystem metabolism by decoupling ecosystem respiration and GPP through limiting light availability, thereby reducing autochthonous C production, and, conversely, by stimulating processing of organic matter through increased surface area (Glover et al., 2019; Honious et al., 2021). The increased surface area of suspended particles in the water column provides microsite habitats for microorganisms (Liu et al., 2013; Ochs et al., 2010), where bacterial production and enzymatic activity are concentrated, contributing substantially to material processing in the water column, particularly in rivers of the fifth order and higher (Gardner and Doyle, 2018; Reisinger et al., 2015). Nutrient dynamics, particularly N, also influence ecosystem respiration, where elevated N concentrations have been linked to increased ecosystem respiration across stream orders (Benstead et al., 2009; Reisinger et al., 2016, 2021; Rosemond et al., 2015). Nitrogen is a key nutrient for microbial growth and is often a limiting nutrient in freshwater rivers (Carroll, 2022). Consistent with this, we found the fastest ER_wc at an agriculturally influenced stream with the greatest TDN and NO${}_{\mathrm{3}}^{-}$ concentrations. Elevated nutrient levels at this site likely stimulate microbial respiration, similar to Cross et al. (2022), who found an increase in heterotrophic respiration in response to N enrichment. Moreover, respiratory processes are typically faster at higher temperatures (Pietikäinen et al., 2005), which can shift riverine ecosystems toward heterotrophy (Song et al., 2018). By stimulating microbial respiration, higher temperatures can also amplify the effects of increasing nutrients (Cross et al., 2022). Collectively, we infer that increasing temperature and nutrients, potentially from anthropogenic inputs, are the most likely drivers of ER_wc in the Yakima River basin. Ultimately, our results emphasize the complex and dynamic roles of the physical, chemical, and biological factors that influence ER_wc in the Yakima River basin and other similar freshwater ecosystems.

4 Conclusions, limitations, and next steps

Our study shows that the ER_wc rates observed in rivers and streams across the Yakima River basin span published rates from studies conducted in rivers across the CONUS and the Amazon River basin. While this study did not measure ER_tot, the observed overlap between ER_wc and literature ER_tot shows the potential relevance of ER_wc to overall stream metabolism. We pose that the high variability observed in ER_wc rates across the basin will likely translate into variable contributions of ER_wc to ER_tot, ranging from negligible to potentially dominant. We anticipate that these influences will not vary systematically moving down the stream network, as we observed a relatively weak association between ER_wc and drainage area across the Yakima River basin. Our results point to more localized control, and the LASSO regression specifically indicated that ER_wc is faster with increasing TDN, stream temperature, DOC, and TSS, consistent with previous work. Overall, our findings show that the complex interactions between physical and chemical factors affect the spatial variability in ER_wc across the Yakima River basin. We encourage future work to expand on our current study by collecting both ER_wc and ER_tot measurements at the basin scale and to consider areal rates to parse the contributions from both the water column and sediments to total ecosystem metabolism.