This study was conducted at two drainage locations on the White Clay Creek, Pennsylvania USA (Figure 1). The upstream site was adjacent to the USGS gage station near London Grove (U.S. Geological Survey #01478100; Coordinates: 39.858 N, 75.783 W). Similarly, the downstream location is next to the USGS gage station #01478245 at Strikersville (Coordinates: 39.747 N, 75.770 W). At the upstream site, White Clay Creek is a second order stream with a drainage area of 750 ha of which 22% is cultivated crops, 43% pasture, 17% forested lands, and 7% is open space developed area (StreamStats/US Geological Survey, 2023). At the downstream location, the creek drains a fourth order watershed with a drainage area of 15,330 ha consisting of 15% cultivated crops, 26% pasture, 24% forested lands, and 30% developed area (StreamStats/US Geological Survey, 2023). Both drainage areas are located within the Piedmont region and have soil types of Glenilg channery silt loam, Worsham silt loam, Wehadkee silt loam, Guthrie silt loam, with Kaolinite clay (Sweeney, 1993), and are underlain mostly by weathered crystalline rocks, shales, schist, gneiss marbles (Pennsylvania Geological Survey, 2017). The potential anthropogenic sources of nitrogen in the watersheds include synthetic fertilizers and manures from agricultural fields and urban areas. The average values for daily discharge are 0.156 (min: 0.066, max: 0.245) m³/s for upstream site and 3.228 (min: 0.651, max: 21.209) m³/s for downstream site (USGS Water Data for the Nation, 2023a,b). The average suspended sediment concentration obtained from our sampling combined with USGS sampled data are 256.4 (min:1.5, max: 1,220) mg/l for the upstream and 1,499.8 (min: 1.6, max: 14,476.8) mg/l the downstream sites.

We collected data on stream water chemistry and N process rates (U_assim and U_den) from selected storms at two drainage locations between September 2020 and December 2021. A total of 33 stream water samples were analyzed for water chemistry and N process rates across the two study sites (20 samples for one April and four September-October events at the downstream site, and 13 samples were used for three September-October events at the upstream site). We used a combination of passive sampling and grab sampling methods to collect stream water on the rising and falling limb of storm events. We built passive samplers following Diehl's design for a modified siphon sampler tower (Diehl, 2008) to collect water on the rising limb. Each sampler tower was built to contain three vertically stacked 1-L HDPE sample bottles inside half of a PVC pipe. Flexible vinyl intake and air tubing were installed through the bottle caps to create a siphon allowing water and SS to flow into the sample bottles once the water level reached the intake tube. We installed metal T-posts to hold the passive samplers near the streambank at each site, which remained at the same location for the duration of the study period. Before storm events, we deployed three passive sampler towers equaling nine total sample bottles at each site with intake tubing facing the stream channel at the same height so that three replicate samples would be collected at a given stage height. At the time of passive sampler deployment at each site, stage height and the distance from the water level to each intake tube were also recorded. This made it possible to calculate the stage height when each sample bottle was filled and determine the sample collection time and discharge values using USGS gage data from each site. On the falling limb of the storm, we removed the passive samplers from the T-posts and collected three replicate grab samples from the water column at each site using a 1-L HDPE bottle attached to a swing sampler (e.g., Nasco 3.65 m Swing Sampler). All samples were then transported to lab facilities in a chilled cooler.

To process samples for collecting water chemistry data, we pooled replicate samples at each stage height from the rising limb and falling limb. We collected subsamples for dissolved organic carbon (DOC) and nitrate-N (NO₃-N) and vacuum-filtered the subsamples through 0.45-μm and 0.7-μm borosilicate fiber filters, respectively. Subsamples for DOC were then analyzed using an Aurora 1,030 W TOC Analyzer (Oceanographic Int., College Station, Texas, United States) and chemical oxidation (Menzel and Vaccaro, 1964). Analysis of NO3−-N was performed via discrete colorimetric analysis using an AQ300 discrete analyzer (SEAL Analytical, Wisconsin, United States) following standard procedures (APHA, 2017). We collected an additional subsample of each pooled sample and vacuum-filtered this through a pre-weighed 0.7-μm borosilicate fiber filter, repeating this step to generate two filters per pooled sample. One set of filter samples we oven-dried at 60°C for 72 hours, weighed on a Sartorius (Goettingen, Germany) MC1 analytical balance, combusted at 500°C for 5 to 6 h (Steinman et al., 2006), and reweighed for calculation of dry mass and ash-free dry mass (e.g., OM). We repeated only the first two drying and weighing steps for the second set of filters and then collected a one-cm subsample from each filter which we weighed and packaged in an aluminum tin capsule. For these subsamples, organic C and N content on SS (as a percentage of dry mass) was determined via combustion using an Elemental Analyzer (University of Maryland, Center for Environmental Science).

We measured ambient denitrification and assimilatory uptake rates for stream water using microcosm experiments. We took three 250-ml aliquots from each pooled sample, generating three replicates per stage height. We poured each aliquot into a microcosm chamber (Kimble 250 ml wide-mouth media bottle) and then enriched all chambers with 1 ml of a 796 mg/l solution of Na¹⁵NO₃. We placed a magnetic stir bar in each chamber, and closed chambers with a screw cap in which septa had been installed. We placed all chambers on a magnetic stir plate and then evacuated (3 min) and flushed with He gas (one min) by inserting tubing with syringe needles attached into the septa of each microcosm chamber (Dodds et al., 2017). The evacuating and flushing cycle was repeated three times before incubation began. We then moved all microcosms to magnetic stir plates inside an incubator set to 25°C. We set stir plates to 360 rpm to ensure particles remained in suspension during the incubation period (Jia et al., 2016). Then, we incubated the microcosms for a period of 24 h in the dark. We collected gas samples at 4 and 24 h after the start of incubation using a gas-tight syringe (Hamilton 25 ml Model 1025TLL) to sample 12 ml of gas from each chamber into 12 ml pre-evacuated Exetainers (Labco Ltd., High Wycombe, United Kingdom). Each gas sample was analyzed for δ¹⁵N of N₂ and N₂O via continuous flow isotope ratio mass spectrometry (IRMS, ThermoScientific Delta V Plus) at the University of California-Davis Stable Isotope Facility. Denitrification rates on SS were measured from a production rate of ²⁹N₂ and ³⁰N₂ following Nielsen (1992). Similarly, assimilatory uptake rates were determined by measuring the increase in the ¹⁵N to ¹⁴N ratio of suspended organic matter at the end of incubation and the tracer 15N: 14N ratio in the microcosm NO₃- following Mulholland et al. (2000). Further details of these denitrification and assimilatory uptake calculations are described in Bacmeister et al. (2022).

We obtained additional information on manual and high frequency (15-min) discharge (Q), Turbidity (Td), and SS concentrations from the USGS stations at the selected drainage locations for the study time period of September 2020 to December 2021 (USGS Water Data for the Nation, 2023a,b). We used this information to develop regression models and compute event mass fluxes for water, SS, and N described below. Additionally, we obtained seven-day average antecedent temperatures of water prior to sampled events at both upstream and downstream sites from the USGS (USGS Water Data for the Nation, 2023a,b) at respective USGS gage stations, and the closest precipitation and 7-day antecedent rainfall data was obtained from “Ag farm weather station” at Newark (National Weather Service, 2023a).

We combined our manually sampled SS data with USGS discharge, turbidity, and SS concentrations to increase sample size and develop stronger regression models. The combined data set included 106 data points for the downstream site (67 on falling limb, 39 on rising limb) and 31 data for the upstream site (25 on falling limb, six on rising limb). We separated rising and falling limbs by the discharge peak, and the end of a storm was identified as the time point when recession discharge was within 10% of the starting storm event discharge (Inamdar et al., 2015). The start of the storm was identified as the point where streamflow discharge exceeded by 10% following the start of rainfall (Inamdar et al., 2015). Separate regression models were developed using a Python programming language for each site and for rising and falling limbs. Rising limb and falling limb datasets were standardized using the natural logarithmic function for developing and testing different statistical models to predict SS concentrations with better r-squared values (R²). We applied simple regression models that are preferable for a relatively small dataset (Trevor et al., 2008; Montgomery et al., 2012). Equations (1)–(4) represent the regression models developed for rising and falling limbs of upstream and downstream sites respectively.

Where SS is suspended sediment concentration (mg/l), Q is discharge (cfs), and Td is turbidity (FNU). Subscripts rv, fv, rd, and fd represent rising limbs and falling limbs for the downstream and upstream sites, respectively.

Three summer/fall storms (Sep-Oct) were available for the upstream site and five storms were sampled at the downstream site across a broader season (Apr-Oct, Table 1). Hydro-meteorological parameters were determined to characterize the events (Table 1). Seven-day average antecedent temperature was determined by determining the average of the daily mean water temperature at the USGS gage stations, and 7-day cumulative antecedent precipitation and event day precipitation were obtained by summing the precipitation data for these time periods. Similarly, we computed the percentage exceedance of peak discharge of storm events based on 2 years of the daily maximum flows (Acharya and Joshi, 2020) at the study sites from 2021 through 2022. Storm discharge for each site was normalized to the drainage area to allow for comparative analysis between the two study sites.

Fluxes for SS were computed by multiplying the flow-weighted average modeled concentration for SS on the rising and falling limbs with the discharge amounts for these periods. The total event SS flux was the sum of the SS fluxes for the rising and falling limbs. We scaled up the U_assim and U_den rates measured in the microcosm experiments to the watershed-level (g/ha) by multiplying microcosm-level rates (g N/kg of sediment/hr) by the duration (hr) and SS flux (kg/ha) of each period (rising and falling limbs) and storm event. Similarly, fluxes for sediment-bound N and C (Sed-N, and Sed-C) were computed in g/ha by multiplying the weighted average %N and %C with the corresponding sediment loads. Other fluxes such as NO₃-N and DOC were calculated based on their weighted average concentration (obtained from grab and passive sampling) and volume of stormflow on rising and falling limbs. The N uptake computed from both study sites were then compared and analyzed against each other, and other similar studies.

Three storm events, E3 (Sep 1), E4 (Sep 23), and E5 (Oct 29) (Table 1), were sampled at both upstream and downstream sites including E3 which is associated with tropical storm Ida (National Weather Service, 2023b). Among the three storms analyzed for the upstream site, storm E4 contributed the smallest runoff (4.1 mm) with a peak discharge of 2.6 m³/s while tropical storm Ida produced a storm runoff of 29.6 mm with a peak discharge of 11 m³/s. The exceedance levels for peak discharges for these events ranged between 1.48–0.28%. For the downstream site, storm runoff, and peak discharge were lowest for E2 (2.5 mm and 22.7 m³/s) and highest for Ida (E3: 40 mm and 203 m³/s). For both sites, antecedent 7-days average stream water temperatures for the Sep-Oct events ranged from 14.5 to 25.4 C and antecedent 7-days accumulated precipitation ranged from 0.80 to 40.1 mm. The spring/April event (E1) had the lowest water temperature (13.8°C) and no antecedent precipitation. Storm flows lasted from 38.8 h (E4) to 72.0 h (E5) at the upstream site, and from 15.5 h (E2) to 65.8 h (E3) for the downstream site. Not surprisingly, tropical storm Ida yielded the largest sediment export (944.1 kg/ha for the downstream site and 215.4 kg/ha for the upstream site) compared to other events (Table 1). This difference is also indicated by the average SS concentrations during storm events—for Ida, they were 525 mg/l and 124.8 mg/l for downstream and upstream sites, respectively, while for other events they ranged from 165.4 to 300.2 mg/l (downstream) and from 44.0 to 71.5 mg/l (upstream).

Carbon exports (DOC and Sed-C) were greater than corresponding nitrogen exports (NO₃-N and Sed-N) for all Sep-Oct storms, however, NO₃-N was slightly higher than DOC for E1. Nutrient fluxes increased with storm size and stream order with tropical storm Ida (E3) yielding the highest fluxes (Table 1). Flow-weighted stream water DOC concentrations were highest for the August-Oct storms, but NO₃-N concentrations were highest for the April (spring) event at the downstream site.

Fluxes of U_assim, and Sed-N increased with size of storm events for both the upstream and downstream sites (Figure 2, Table 1). The increase in these fluxes, were however, not directly proportional to the magnitude of the events. Flux values for U_den were more variable than U_assim with respect to storm event magnitude, especially for the downstream site (Figure 2). Importantly, storm event Sed-N exports were 1–2 orders of magnitude greater than the corresponding U_den and U_assim fluxes (Table 1). As a % of the Sed-N exports, U_assim fluxes were consistently less than U_den for the upstream site and ranged from 4.1–22.4%. A mixed response was observed for the downstream sites with U_assim fluxes ranging between 0.6–8.8%. U_den varied between 0.1–40.6% for the study sites (Table 1). The April event (E1) measured at the downstream site yielded the lowest process fluxes.

Storm event Sed-N exports at the downstream site were 2–4 times the values measured at the upstream location (Figure 3). In comparison, the responses for N process fluxes (U_den and U_assim) were mixed (Table 1, Figure 3). Assimilatory N (U_assim) fluxes were greater at the downstream vs. the upstream site for events E3 (three times greater) and E4 (1.4 times greater), but similar for E5. In contrast, U_den flux was greater at the upstream site for two (E3 and E5) of the three common events. Tropical storm Ida in particular, produced a denitrification flux that was four times greater at the upstream vs. the downstream site (Table 1, Figure 3). U_den flux for the downstream site for E5 was unexpectedly very low (close to zero; Table 1), with low U_den rates in the rising limb and non-detectable denitrification during the falling limb of the storm.

N process fluxes were consistently higher on the hydrograph falling limbs vs. the rising limbs regardless of storm size and drainage location (Figure 4). While Sed-N exports were also higher on falling limbs except for the tropical storm Ida at the upstream site, the differences between rising and falling for U_den and U_assim are much greater than for Sed-N. At the downstream site, U_den, and U_assim fluxes before peak flow were on average 71.4–99.8% and 98.6–99.8% smaller than those observed during the falling limb, respectively (Figure 4, Upstream). Similarly, rising limb values for U_den and U_assim fluxes at the upstream site were on average 95.3–99.8% and 96.9–99.9% smaller than mean fluxes in the falling limb (Figure 4, Downstream).

We found that while the storm-event Sed-N exports were consistently greater at the downstream site, the same was not true for U_den and U_assim fluxes. U_assim fluxes were higher for downstream site for two events and similar for the third one. However, U_den (for two of the three events) were greater for the upstream drainage location. Previous studies have typically reported the rates, rather than watershed-scale fluxes [in M L²] for denitrification (Christensen and Sørensen, 1988; Royer et al., 2004; Arango et al., 2008; Reisinger et al., 2015; Xia et al., 2017b; Wang et al., 2022). To allow for comparative analysis and evaluation with previous studies, we computed areal rates of U_den for our sites (Table 2). This was done by multiplying our volumetric rates with the average high flow depth estimated empirically from measured depths from USGS gage stations at those sites (we estimated high flow depth for the upstream site to be ~1.3 m and ~1.15 m for the downstream site). Our U_assim rates (0 to 5.82 mg m⁻³ hr⁻¹ for the downstream site, and 0 to 2.32 mg m⁻³ hr⁻¹ for the upstream site) were within the range of rates (0.001–363 mg m⁻³ hr⁻¹) reported by Reisinger et al. (2015) across first to fifth order watersheds with contrasting landscapes in mid-West US.

Our water column U_den rates for both sites (0.00 to 6.08 mg m⁻³ hr⁻¹ for the downstream site and 0.00 to 9.82 mg m⁻³ hr⁻¹ for the upstream site) were within the broad range (−3.20 to 22.80 mg m⁻³ hr⁻¹) of those reported by Wang et al. (2022) for third to eighth order rivers in China. Wang et al. (2022) found that denitrification fluxes increased with stream order or drainage area and attributed it to increase in suspended sediment, total organic carbon, and water column depth. Reduced dissolved oxygen and increased microbial communities associated with the suspended sediments have also been identified as key factors in enhancing water column denitrification removal (Jia et al., 2016; Xia et al., 2018; Wang et al., 2022). In another study in China on the Yellow River, Xia et al. (2021) found that N₂ and N₂O fluxes from the river were highest in the middle reaches but lower in the upper and lower reaches (specific stream orders corresponding to “upper” and “lower” reaches were not provided). Xia et al. (2021) attributed the elevated N₂ and N₂O fluxes in the middle reach to the increased suspended sediment concentrations in this reach. Reisinger et al. (2016) also reported variable denitrification rates (range: 0.00 to 4.90 mg m⁻³ hr⁻¹) among five mid-Western rivers in the US.

Our denitrification rates were closer to the rates obtained by Reisinger et al. (2016) and did not increase with an increase in stream order (the upstream site had higher U_den rates compared to the downstream site). This was surprising given that suspended sediment and DOC concentrations were greater at the downstream vs. the upstream site while dissolved NO₃-N concentrations were comparable (Table 1). We speculate that greater fraction of agricultural land use providing more labile organic carbon and nitrogen in the upper watershed coupled with finer sediments could have played a potential role in the elevated denitrification rates/fluxes measured for the upstream site. Grain size distribution, was however, not available in our study. We also recognize here that these drainage scale comparisons were based on limited data from only three storm events (all in Sep-Oct time frame; Figure 3) and additional events across multiple years and multiple drainage locations are needed to make a robust assessment for the influence of drainage size on denitrification fluxes. Contrary to U_den, U_assim fluxes were however higher at the downstream site (Figure 3) and might be due to factors such as increased light, temperature, and algal biomass in the wider, downstream reaches (Rode et al., 2016). Indeed, mean stream water temperature was consistently higher at the downstream vs. the upstream site (Table 1).

Sed-N and assimilatory uptake fluxes displayed a pronounced increase with storm event magnitude (Figure 2), although the increase in Sed-N was greater than U_assim. The sharp increases in sediment bound N were not surprising given that the mobilization of suspended sediments and associated particulate N in storms is well recognized (Royer et al., 2006; Edwards and Withers, 2008; Inamdar et al., 2015; Jiang et al., 2020). The elevated values of U_assim and U_den during storms suggest that these fluxes are not trivial and need to be accounted for in reach and watershed scale N budgets and assessments. What was surprising though was that while U_den increased with storm magnitude at the upstream site, the same pattern was not repeated for the downstream site (Figure 2). Previous studies have typically observed U_den to increase with runoff amounts/depths and suspended sediment concentrations that provide valuable sediment surface area for microbial habitat and reactions that may facilitate denitrification and assimilatory processes (Xia et al., 2017b; Wang et al., 2022). Increased suspended sediment concentrations have also been found to enhance anoxic microsites that may support increasing denitrification loss (Beaulieu et al., 2009; Yu et al., 2013; Xia et al., 2021). Thus, the variable storm response at the downstream site suggests that additional unknown factors could influence the U_den process in the water column.

Temperature is also an important control that could influence the microbial growth processes responsible for water column assimilatory and denitrification fluxes (Xia et al., 2021; Wang et al., 2022). Elevated temperatures could also decrease the amount of dissolved oxygen in the water column making the environment more suitable for low oxygen nitrate reduction processes. Xia et al. (2021) reported higher N₂ and N₂O fluxes for summer vs. spring for the river network. While we did not specifically focus on seasonality and temperature differences across storms, our observations for the downstream site (Table 1) indicated that the U_den fluxes for the April/spring (E1) event were lower than those measured for most of the summer storms (except E5). Our sample size is however very small and additional storms across seasons need to be measured for a more thorough assessment.

This study provided important new insights on within-event patterns of (rising vs. falling hydrograph limb) U_den and U_assim fluxes which have not been reported before. Both N fluxes were greater on the falling limb of the discharge hydrograph as opposed to the rising limb (Figure 4). Duration and runoff volumes were both greater during the falling limbs of the storm hydrographs (Table 1) and likely contributed to the larger fluxes for this period. Previous studies have shown that hydrologic flow paths, sediment sources and characteristics, and nutrient concentrations could differ on the rising vs. falling limbs of the hydrograph (Buttle, 1994; McGlynn and McDonnell, 2003; Inamdar et al., 2013; Klaus and Mcdonnell, 2013; Johnson et al., 2018). Rainfall input and high erosive energy of runoff on the rising limb typically sustain greater grain size for suspended sediments vs. the falling limb (Walling et al., 2000). The flow turbulence on the rising limb could also facilitate well-mixed and more oxygenated conditions in the water column vs. the less turbulent flow on the recession limb of the hydrograph (Tabarestani and Zarrati, 2015). DOC, which is important energy source for heterotrophic denitrification, has also been typically reported to peak following the discharge peak (Inamdar et al., 2011; Dhillon and Inamdar, 2014). In a previous study in our upstream watershed in White Clay Creek (Kan, 2018), microbial communities involved in N processes were observed to vary across the storm hydrograph. Kan (2018) reported different microbial compositions as well as greater abundances of nitrification and denitrification genes on the recession limb vs. the rising limb of the discharge hydrograph. Terrestrial microbes flushing into a stream during storm flow might need time to adapt and thrive, showing up on the falling limb. We hypothesize that all or some of these factors likely contributed to the higher amounts of N fluxes that we observed on the falling limb of the discharge hydrograph. Future studies are needed to explicitly determine the influences of these factors on the N fluxes.

Our results show that while the sediment bound N exports were much greater than U_assim and U_den during the large storms, the amounts of U_assim and U_den were not trivial and watershed N budgets and models would benefit from incorporating water column removal processes. Currently, most watershed-scale models for nutrient budgeting and management, including the more popular ones like the Soil Water Assessment Tool (SWAT; Arnold et al., 2012), and Storm Water Management Model (SWMM; Rossman and Huber, 2015) do not include water column denitrification processes (assimilatory uptake is simulated to various extent; Arnold et al., 2012; Rossman and Huber, 2015). This is an important knowledge gap. Models like SWAT and SWMM incorporate erosion and sediment transport algorithms that can describe sediment and sediment bound nutrient exports from the land surface into the stream (Arnold et al., 2012). The SWAT model already characterizes (via the In-channel module)—stream/river water depth and travel time for simulated reaches, dissolved oxygen, and suspended sediment concentration, transport, and particle size distribution. The model also has limited algorithms to characterize bacterial biomass (but not type) associated with channel sediments. These algorithms and sediment parameters can be leveraged and simplified equations linking water column denitrification to parameters such as water temperature, dissolved inorganic N, organic C, and suspended concentrations can be included. Such equations based on empirical data have already been developed by Xia et al. (2017b), Pang et al. (2022), and Wang et al. (2022) for various rivers in China. However, regression parameters for these equations likely vary with physiographic regions and climate and so region/climate specific values may be needed. For example, a comparison of our measured values (and fitted regression equations) against predicted estimates based on equations developed by Xia et al. (2017b) indicate that our measured fluxes were under predicted (Figure 5). Thus, accurate estimates of water column denitrification may be needed if we want to develop a robust estimate of these processes at the reach and watershed scales.

We also recognize that N process rates may vary with reach and watershed scales and algorithms and equations developed at the small mesocosm/microcosm scales may not necessarily apply as is at the larger scales. Addressing scaling issues for processes and their parameters has been an ongoing challenge across multiple disciplines (e.g., Blöschl and Sivapalan, 1995). This challenge stems from the inability to measure process rates and parameters at large scales. This applies here too given the practical challenges with measuring U_den and U_assim at the scale of the full stream reach. Currently, the best strategy in such situations is to collect as many N process measurements as possible and for a wide range of conditions (including discharges, suspended sediment concentrations, water temperatures, and dissolved oxygen) and assume that the relationships provide some representation of processes at the larger scales.

Additionally, climate change could further increase the frequency and intensity of the largest storms (Wuebbles et al., 2017). Increased intensity of storms is expected to increase sediments and sediment-bound nutrients from watersheds (Dhillon and Inamdar, 2013). This could increase the potential for water column denitrification losses as well as assimilatory uptake, particularly in combination with increasing air and water temperatures, thereby adding additional urgency to improving N fluxes in watershed budgets and models.