The experiments were carried out in a 353-ha tidal tributary to the Chesapeake Bay found in Anne Arundel County, Maryland (Figure 2). Rock Creek’s watershed is 80% residential and 20% forested. Due to the poor water quality over the past few decades (e.g., anoxia and extensive algal blooms), the county installed an engineered destratification system in 1988 to bring dissolved oxygen back to former levels (Harris et al., 2015). Herein, we call this the “large bubble aeration” system. Up until 2019, this system was made up of 138 ultra-coarse air diffusers distributed along a pipe that lines the middle creek channel with ∼20 mm bubble size (CH2M_Hill, 2011), as described in design specifications reported by Dames and Moore (1988). The goal of the system was to vigorously overturn the water column in order to continuously introduce oxygen via de-stratification into the bottom waters, and the system was run continuously throughout the day. Every year, the aerators are turned on June 1 and turned off October 1, in order to minimize the effects of summertime hypoxia. An early study of this system determined that the zone of aeration influence was ∼74 ha and that bottom waters remained oxic when diffusers were ON, but became anoxic within one tidal cycle when aerators were OFF (Harris et al., 2015; Harris et al., 2016). In spring 2019, the aeration system was upgraded to fine bubble diffusers to provide more oxygen to the water column, herein referred to as the “small bubble aeration” system. The goal of the diffusers is primarily to overturn the water to allow re-aeration to occur at the water surface, not to add oxygen from the bubbles themselves. There are two 213 m long diffusers emanating from the shore mainline, which provide air at a rate of 180 scfm (standard cubic feet per minute) in a continuous bubble pattern (0.26 scfm ft⁻¹) of bubbles 3 mm in diameter. The surface water expression of the new aeration system is shown in Figure 2B. With this new system, the county public works turns OFF the aerators every night to reduce neighborhood noise and energy consumption.

Over the course of 4 years (2016, 2018, 2019, 2021), the water column and sediments were sampled along the creek both within and outside the aeration zone and with both engineering designs (Tables 1, 2). Stations within the aeration zone included RC1, which is located at the up-creek limit of aeration, and RC2 which is found mid-channel and directly in the aeration zone (Figure 2). RC7 is ∼1 km downstream from the end of aeration, still within the zone of influence of aeration and where the creek widens, and RC9b is a background site, close to the mouth of the Patapsco River and outside the zone of influence of aeration (Figure 2). These stations were introduced in previous studies (Harris et al., 2015). Water depths are between 1.5–3.5 m. In 2021, four upstream stations were added to determine the river influence to the aeration zone (Table 2).

To study the effects of the aeration, our experimental approach was to sample while the aerators had been on for about 1 month, referred to as the “ON” sampling event, which occurred during the daylight hours. In the evening of the “ON” sampling, the aerators were turned off and waters sampled 1 day later (2016; large bubble aeration), 7 days later (2018; large bubble aeration), 13 days later (2019; small bubble aeration), and 3 days later (2021; small bubble aeration); these are referred to as the “OFF” sampling events (Table 2). After completion of the experiments, the aerators were turned back ON for the remainder of the season. By sampling at different time periods after aerators were turned OFF, the experiment addressed the question of the impact aeration had on methane flux from the Rock Creek estuary.

During each field campaign, we collected water column hydrographic and chemical profiles, discrete water samples from the surface and bottom depths, and shallow (∼30 cm) sediment cores via small boat. Water column temperature, salinity, and dissolved oxygen levels were recorded with a YSI EXO2 multiparameter sonde. In 2018, a benthic lander was placed at a dock (location shown as black cross in Figure 2) at the edge of the aeration zone that included continuous temperature, salinity, and dissolved oxygen sensors (YSI EXO2). Salinity is reported in practical salinity scale which has no units. Wind speeds were determined from a handheld anemometer (Weatherhawk, Windmate WM-200).

For water column and sediment vial samples, the headspace was equilibrated with the dissolved methane from the aqueous sample and the headspace analyzed for the ppmv methane. For all samples, an aliquot of the headspace was extracted from the vials and injected onto a gas chromatograph (SRI 8610C multi-gas) equipped with a HayeSep D packed column and a Flame Ionization Detector in order to quantify methane concentrations (Magen et al., 2014). Certified standards (Airgas, Inc.) were used for the calibration curve. Analytical precision is 3% and all measurements were above detection limit of 2 ppmv. Resultant partial pressures were then used to calculate dissolved methane concentrations (in nM) in either the water column or the porewater using Henry’s law according to equations in Magen et al. (2014) and porosity corrections according to Lapham et al. (2008a).

For sulfate and chloride concentrations in the sediment porewater samples, the tube containing whole sediment was centrifuged (3000 RPM, 30 min, 20°C, Sorvall^© RT 6000D) and the resultant supernatant filtered with a 0.2 µm syringe filter. Samples were then diluted (1:135) in Milli-Q water and analyzed on a Dionex ICS 1000 ion chromatograph (IonPac AG22 4 × 50 mm guard column, IonPac AS22 4 × 250 mm analytical column, and ASRS 300 4 mm suppressor) with an AS40 Autosampler. Water samples from the 50 cm OsmoSampler sections were also measured for chloride with the same dilution to calculate salinity. Certified IAPSO seawater standard (Ocean Scientific International Ltd.) was used for the calibration curve. Analytical precision is 2% and all measurements were above detection limit of 0.05 mM for sulfate.

Air syringes were directly connected to the intake of a cavity ring down spectrometer (CRDS, Picarro 2201i) to measure for methane concentrations and methane stable carbon isotopes. For the water samples, 10 mL of degassed brine was added to the vials to displace the headspace and injected into the small sample isotope module (Picarro, Inc.) to introduce a small sample to the CRDS, similar to procedure in McIntosh Marcek et al. (2021). Isotope values were obtained through calibration with three Vienna Pee Dee Belemnite (VPDB) referenced standards (−23.9‰, −38.3‰, and −66.5‰ (±0.2‰); Isometric Instruments). Isotopic results are reported using the δ¹³C notation in per mil (‰), where δ¹³C = (R_sample/R_standard -1)*1,000 and R = ¹³C/¹²C. Analytical precision is 2% for concentrations and 4‰ for stable carbon isotope ratios.

The air-water flux of CH₄, F, was determined for all discrete sampling campaigns using the updated flux equations presented Wanninkhof (2014): F = k ( C w − C e q ) (1) where k is the gas transfer velocity (length time⁻¹), C_w is the measured surface water concentration, and C_eq is the CH₄ concentration in equilibrium with the atmosphere at in situ conditions (Yamamoto et al., 1976). There are several versions of Eq. 1, mostly based on wind speed. Here we employ the formulation and parameterization of Myllykangas et al. (2020), which reports a k value adapted from Raymond and Cole (2001): k = 1.91   e 0.35 u ( S c 600 ) − 0.5 (2) where u is the average wind speed and Sc is the Schmidt number for CH₄ in freshwater calculated from Wanninkhof (2014). Wind speeds were obtained from a nearby NOAA buoy (National Data Buoy Center BLTM2, Baltimore, MD) and averaged over the 3 days prior to sampling. Since the buoy is located 15 km away from Rock Creek, we compared the buoy to the handheld anemometer readings and found they compared within 7%. The buoy wind speed was used for all stations. Wind speeds varied between 2–3.2 m s⁻¹, which translated to k values varying between 3.8 and 5.9 cm h⁻¹, similar to values found in MacIntyre et al. (2010) in a lake system and mangrove dominated estuaries (Rosentreter et al., 2017). Air-water fluxes were then calculated using Eq. 1 and reported as µmol CH₄ m⁻² d⁻¹.

Since the calculated air-water flux is inherently based on assumptions of a stagnant boundary layer, the calculated air-water methane flux when the aerators are ON will be underestimated. To constrain this better when the aerators were ON, we calculated an air-water methane flux from direct bubble transport to surface water by applying the bubble radius of the system (3 mm, Mobley Engineering, Inc., personal communication) to an existing bubble model output to determine the mass transfer coefficient (Figure 4 in McGinnis and Little, 2002). This mass transfer coefficient for the bubble radius in the system is 0.04 cm s⁻¹ (or 144 cm h⁻¹) which was then used in Eq. 1 to calculate a modified air-water methane flux for 2018 and 2019 at RC1 and RC2.

The sediment-water methane diffusive flux was calculated from Fick’s first law: J C H 4 − S W I =   − φ D s d C d x (3) where J_CH4-SWI is the methane flux (µmol CH₄ cm⁻² yr⁻¹) at the sediment-water interface, φ is the porosity (0.8), D_s is the sedimentary methane diffusion coefficient (cm² s⁻¹), x is the vertical sediment depth (cm) and dC/dx is the concentration gradient of methane. D_s was calculated for each station, corrected for tortuosity and in situ pressures (based on water depth), temperatures, and salinity (Millero, 1996), and was ∼1.7 × 10^–6 cm² s⁻¹. The gradient term was calculated between the uppermost porewater measurements; which usually started at 1.5 cm into the sediment. Thus, this gradient is most likely overestimated because it ignores any oxidation processes that might occur in that upper 1.5 cm of sediment. Using this gradient, the diffusive fluxes were then calculated for each station and time point with Eq. 3. For convenience, fluxes are reported as positive but represent flux out of sediment).

We applied a simple box model to the flux data with the assumption that the sediments are the only source of methane to the water, and the atmospheric flux was the only sink of methane in the water. If the two balanced, then there would be no additional contributions to the methane budget in Rock Creek. If the fluxes did not balance, we could invoke additional microbial oxidation or production in the water column and/or advective transport of methane from up- or down-stream. To do this, at each station for 2018 and 2019, we subtracted the atmospheric methane flux from the sedimentary methane flux, then assigned the net difference as the water column methane inventory anomaly. For the aeration sites (RC1 and RC2), both air-water and sedimentary fluxes were averaged together. We should note that when the aerators are ON, we used the fluxes estimated from the stagnant boundary layer model which most likely underestimates the atmospheric flux.

Student t-test was used in Excel (two-tailed, paired) to determine if the methane concentrations calculated from the discrete sampling events and the OsmoSampler samples were significantly similar. To create the paired methane dataset for this test, we extracted the OsmoSampler methane concentration that was sampled at the same time as the discrete methane measurement from the water column at each station. Please note that the OsmoSamplers average over 6–12 h. This comparison was done with the fast OsmoSamplers during 2018 (see Section 3.6 for result).

The experimental design was intended to measure methane concentrations in well oxygenated, bubble-influenced waters when aerators were ON, and then hypoxic or anoxic waters without bubble transport when the aerators were turned OFF. Based on previous experiments in Rock Creek, the change to hypoxic conditions occurred within 1 day of turning off the aerators (Harris et al., 2015). However, for our manipulations, the goal was to observe changes over the longer term (up to 13 days) and the systems response. It is important to remember there was a 7-day difference between ON and OFF treatments in 2018 and 13-day difference in 2019.

In 2018 during aeration, the dissolved oxygen (DO) concentrations were ∼7 mg L⁻¹ at RC1 and RC2, and well mixed (Figure 3 top panel). At RC7, the water column was still oxygenated (>4 mg L⁻¹; comparable to conditions at RC9b). After aerators were OFF for 7 days, RC1 and RC2 became hypoxic throughout the water column below 0.5 m, and the other stations had weakly stratified water columns, but with hypoxic bottom waters below 2.5 m. Dissolved inorganic nitrogen concentrations averaged 14.4 µM (ranging between 0.8 and 19 µM for all stations) when aerators were ON and 7.2 µM (ranging between 3.7 and 12 µM) when aerators were OFF (Table 2). Salinity was ∼5–5.5 (Figure 3). Water temperatures in 2018 varied between 26–30°C, and were similar for both ON and OFF treatments (data not shown).

In 2019, dissolved oxygen was relatively high in surface waters when the aerators were OFF at RC1 and RC2, giving way to oxygen-depleted conditions below 2 m (Figure 3). Dissolved inorganic nitrogen concentrations averaged 27.8 µM (ranging between 5 and 36 µM for all stations) when aerators were ON and 4.9 µM (ranging between 1 and 13 µM) when aerators were OFF (Table 2). Salinity was lower in 2019 than in 2018 and was nearly 3.5 when aerators were ON, and 4.5 when aerators were OFF (Figure 3). Water temperatures were between 26 and 29°C during the ON treatment and were warmer during OFF treatment (data not shown).

The tidal stage at each station varied over the discrete sampling time points (Supplementary Figure S3). In 2018 when aerators were ON, RC7 was sampled first at the ebbing tide, RC9b at low tide, and RC1 and RC2 at a high tide. When aerators were OFF, RC1 and RC2 were collected close to high tide or when waters were just beginning to ebb. RC7 and RC9b were sampled on the ebb tide. In 2019 when aerators were ON, RC9b, RC1, and RC2 were collected on flooding tide, and RC7 was collected right after high tide. When aerators were OFF, all stations were collected near the high tide.

Dissolved methane concentrations in surface and bottom water of Rock Creek varied over space and time (Figure 4). Overall, concentrations ranged between 150 and 1,500 nM, orders of magnitude higher than atmospheric equilibrium (which is ∼3 nM), and were higher at stations RC1 and RC2, within the aeration zone, compared to stations closer to the Patapsco River (RC7, RC9b). In 2016, concentrations at RC2 and RC7 were around 400 nM throughout the period of measurements, regardless of aeration status (Figure 4A). In 2018, during aeration, the waters were relatively well mixed between surface and bottom waters (Figure 4B). Once aerators were turned OFF, there was an increase in bottom water methane at RC1 and RC2, and not much change at RC7 and RC9b. After 7 days, the surface waters were enriched in methane compared to the bottom water. In 2019, the concentrations between surface and bottom water followed expectations: during aeration, the water column was well mixed so there was little difference between surface and bottom waters and when aerators were turned OFF, methane concentrations were higher in the bottom water than surface waters after 13 days (Figure 4C).

The 2021 field campaign was designed to resolve the upstream contribution of methane to the aeration zone. Overall, methane concentrations were lower than 2019 values but also showed the same pattern of higher concentrations at RC1 and RC2, than at RC7 (Figure 4D). At RC7, concentrations were similar between surface and bottom, and regardless of aeration status. At RC1 and RC2, methane concentrations were lower when aerators were OFF, which was unexpected. The other unexpected result was to record higher methane concentrations in surface waters than bottom waters during both ON and OFF periods (Figure 4D). Measurements made upstream of the aerators, which were only made in 2021, showed that surface waters were always higher than the bottom waters and that concentrations during the ON treatment were always higher than the OFF. Furthermore, methane concentrations were highest in the most upstream station, and declined downstream and into the estuary, such that water flowing into the aeration zone from upstream were enriched with methane relative to the aeration zone itself.

The water column methane concentrations are also presented in Figures 5A–D, 6A–D as compilation figures showing the water column and sediments in a holistic view. Here we add to the water column concentration data the stable isotopic ratio of methane carbon to distinguish source of this methane (Figures 5I–L, 6I–L). Regardless of station or aeration status, δ¹³C-CH₄ values ranged from −69 to −51‰, with an average value of −61.4 ± 3.6‰, which is near the standard deviation of the method.

Methane concentrations measured from the sediment porewaters increased with sediment depth (Figures 5, 6, brown colored panels). Surface concentrations were at ∼µM levels and increased to as high as 1,300 µM at the bottom of the core. In both years, concentrations were higher in the aeration zone (RC1 and RC2) and outside the aeration zone (RC7) compared to the background site (RC9b). Yet, there is variability in the sediment profiles. For example, at RC1 and RC2, sediment methane concentrations were lower in 2018 than in 2019, regardless of aeration status. After 13 days of no aeration, methane concentrations were higher in the sediments at RC2 and RC7. As noted, Figures 5, 6 also contain methane water column parameters for comparison purposes.

Methane increases in sediment porewaters are typically associated with a drawdown of sulfate in the surficial depths due to sulfate reducers outcompeting methanogens for substrates (Hoehler et al., 1994). Therefore, we also measured sulfate in porewater to help our understanding of anaerobic biogeochemical processing. Sulfate concentrations decreased downcore in all stations except RC9b although the depth of low sulfate (SO₄ < 0.5 mM) varied (Supplementary Figure S4). The depth of low sulfate typically coincides with the increase of methane, and is known as the sulfate methane transition (SMT) depth. The SMT is an area of active anaerobic methane oxidation via sulfate reduction (Jørgensen et al., 2020) and is a useful metric to show how active the anaerobic microbial community is in a sediment column (Figures 5M–O, 6M–O). While we expected to see the SMT depth shoal when aeration was turned OFF, there was no consistent pattern of the depth of the SMT with aeration status for both years, although we will specifically present SMT depth patterns in each year below. Using the gradients from the top of the cores, the methane flux to the sediment-water interface varied across space and time, and ranged between 0.1 and 700 µmol m⁻² d⁻¹ (Supplementary Figure S5).

Chloride concentrations were also measured as a conservative tracer and as a way to validate any depletion of sulfate coming from sulfate reduction and not a consequence of groundwater. In 2018, while the chloride concentrations showed a slight increase in depth, the depth averages are as follows: 76 ± 8 mM (RC1, RC2), 100 ± 14 mM (RC7), and 85 ± 10 mM (RC9). Since these chloride values are within the range of what would be expected given the overlying water salinity, we represent any conservative mixing in terms of how it might affect sulfate concentrations. Given the rule of constant proportions, we calculated the range of sulfate values that would be estimated given those chloride concentrations (shaded rectangles in Supplementary Figure S4).

Methane stable carbon isotope ratios were measured to help distinguish the fate of methane formed in the sediment. Our measurements revealed two patterns in δ¹³C-CH₄ values with depth: 1) δ¹³C-CH₄ increasing with depth and 2) δ¹³C-CH₄ peaks at intermediate depths associated with the SMT. In 2018, at RC1, the δ¹³C-CH₄ values were similar during the ON and OFF conditions in that the surface was ¹³C depleted (between −80 and −70‰), they became heavier with depth to as high as −50‰, and then decreased again to near surface sediment values (Figure 5M). The depth of the SMT deepened with aerators OFF. At RC2, during the ON condition, δ¹³C-CH₄ values were ∼−80‰ and increased with depth in the core (Figure 5N). During the OFF condition, δ¹³C-CH₄ values were around −60‰ in the shallow depths and quickly decreased to −80‰ at the SMT. At RC7, just outside the aeration zone, δ¹³C-CH₄ values were ∼−70‰ at the surface when waters were aerated and decreased to −85‰ at the bottom of the core (Figure 5O). During the OFF situation, values showed a similar trend at the surface but then showed a mid-depth minimum of −80‰ at the SMT.

In 2019, the δ¹³C-CH₄ value trends showed more consistency across the stations (Figures 6M–O). At RC1, values decreased at the surface from as high as −40 to ∼ −70‰ where values remained for about 15 cm into the sediments (Figure 6M). RC2 and RC7 show almost the same isotope profiles where values decrease downcore, but the δ¹³C-CH₄ values are ∼10–15‰ higher when aerators were ON compared to when they were OFF (Figures 6N,O).

Methane concentrations were measured in the air above each station during 2018 and 2019. In 2018, the average atmospheric methane concentration across all sites was 1.84 ± 0.06 ppmv, and didn’t vary between ON and OFF conditions, except at the dawn sampling (Figures 7A,B). The average δ¹³C-CH₄ value of the background methane was −51.4 ± 10‰ (Figures 7C,D). Dawn sampling on the ON and OFF days showed elevated methane concentrations, reaching as high as 3 ppmv at RC7 which had a δ¹³C-CH₄ value of −90‰ (Figures 7A,B).

To calculate the air-water methane flux, we used two approaches. The first used the stagnant boundary layer model and most likely underestimates the flux for when the aerators are ON. Using this model, the air-water methane fluxes ranged between 300 and 1,500 µmol CH₄ m⁻² d⁻¹ (Figure 8). The flux was higher at RC1 and RC2 than other stations, regardless of aeration status or year. In 2018, the flux at the RC1 and RC2 was higher when aerators were ON after 7 days, whereas in 2019, the flux was lower when the aerators were ON. The second approach was only carried out when the aerators were ON and assumed methane was being stripped from the water as the aerator bubbles traveled up the water column. The calculated fluxes were much higher than the fluxes from the stagnant boundary layer (Figure 8 extended arrows to blue dots). In 2018, at RC1 and RC2, air-water methane flux was 30,730 and 19,380 µmol CH₄ m⁻² d⁻¹, respectively. In 2019, at RC1 and RC2, air-water methane flux was 14,669 and 8,342 µmol CH₄ m⁻² d⁻¹, respectively.

The sediment and air-water fluxes were used in the box model to determine if there are additional sources or sinks of methane beyond what is coming from the sediments and being lost to the atmosphere (Figure 9). In 2018, during large bubble aeration, there was a large source of methane (positive values in Figure 9) at the aerators, and actually a methane sink from RC7 when the aerators were off. In 2019, there was a methane source across all stations, regardless of aeration status, and this source was fairly constant across the sites (Figure 9). The exception to this was at the aerators when they were ON; there was a large methane sink (negative value in Figure 9). However, since the atmospheric flux during the ON status is mostly likely underestimated, this exception is most likely a methane source too.

Using OsmoSamplers, two separate records of dissolved methane concentrations from bottom water were obtained. The first was from the OsmoSampler deployment in 2018 that spanned 1-week using fast pumps with ∼1 day resolution (Figure 10). The time stamps assigned were verified by comparing salinity (as calculated from chloride concentrations) in the OsmoSampler coils and the salinity from the sensor packages (Supplementary Figure S6). The salinity comparison shows relatively good agreement, especially at station RC7. The overall trend is similar between the sensor and OsmoSamplers at RC1 and RC2, but as we have observed in previous studies, the absolute salinity values did not match well at these stations (Gelesh et al., 2016). The highest methane concentrations came 3 days after the aerators were turned OFF at RC2 and reached almost 3,000 nM (Figure 10). Concentrations were also high at RC1 during this time. The timing of this methane peak came right after an event where dissolved oxygen (measured between stations RC2, and RC7) increased rapidly to ∼8 mg L⁻¹ oxygen (Figure 10). After 14 July 2018, methane concentrations decreased to less than 1,000 nM and were similar at all stations. Methane concentrations from OsmoSamplers were cross-checked with our discrete samples and we see no statistical difference between the two (p = 0.72); which is the first time this has been verified in field tests. Methane concentrations at RC7 remained lower than the other stations.

The second time-series record of methane concentration came from OsmoSamplers deployed through the summer and into the fall of 2018 and 2019, and contained slow pumps that give ∼ weekly resolution (Figure 11). The temporal pattern was not the same each year. In 2018, at RC1 and RC2, the initial concentrations before aeration were lower (∼400 nM), and then almost doubled when aerators were turned OFF (Figure 11A). Once they were turned back ON after our experiment, concentrations at RC1 continued to decrease at a rate of ∼13 nM day⁻¹ (linear fit with R ² = 0.8); whereas at RC2, concentrations continued to increase through July and finally peak in August at 17,000 nM. Concentrations at RC2 then decreased and reached ∼1,000 nM for the remainder of the timeseries. At RC7, concentrations didn’t show much change with time and averaged 842 ± 265 nM. In 2019, we captured much higher temporal resolution with the samplers which started about 2 weeks before our experiment began (Figure 11B). Overall, concentrations were lower than in 2018 and ranged between 110 and 1,667 nM. There were concentration differences across sites, where methane concentrations at RC1 averaged 368 ± 100 nM; RC2 averaged 558 ± 136 nM; and RC7 averaged 400 ± 270 nM (Figure 11B). The bottom water temperature varied between 23–28°C in 2018 with some variability (Supplmentary Figure S7), whereas in 2019, the temperature gradually increased from ∼23°C to a high of ∼30°C in August and then decreased into the fall where a sudden decreased to less than 20°C occurred when the aerators were turned off (Supplementary Figure S7).

We hypothesized that the air-water methane flux would be higher during aeration then when the aerators were OFF. When we simply apply the stagnant boundary layer model to calculate the fluxes, we see that sites within the aeration zone (RC1 and RC2) had higher methane fluxes than downstream, regardless of aeration status (Figure 8). Yet these fluxes are most likely underestimates when the aerators are ON, as we see with the modified flux calculation (Figure 8, extended arrows to blue dots). We further hypothesized that under small bubble aeration, the flux would be lower than under large bubble aeration (Figure 1); which is what the fluxes showed in 2018 (large bubble aeration) versus 2019 (fine bubble aeration, Figure 8), supporting this hypothesis. The other observation was that in 2019, when the fine bubble aeration was installed, the air-water methane flux was lower during aeration then when the aerators were turned OFF. We have already stated that this flux is most likely underestimated. Future work would benefit from directly measuring this flux with floating chambers to more precisely quantify this flux.

Rock Creek is a source of atmospheric methane, regardless of aeration status, or site. The air-water methane flux from Rock Creek varied between 0.2 and 1.5 mmol m⁻² d⁻¹ (note change in units to compare to literature values), which was similar to fluxes measured from several estuaries (Table 3), and higher than those from oceanic environments, which vary between 0.0001 and 0.1 mmol m⁻² d⁻¹ (Bižić et al., 2020). Rock Creek methane fluxes are on par with a shallow subarctic lake which reached almost 0.4 mmol m⁻² d⁻¹ (Jansen et al., 2020), even though at times of ice-out, these lakes can release as much as 75 mmol m⁻² d⁻¹ (McIntosh Marcek et al., 2021). Surface water concentrations were also similar to those measured from an aerated eutrophic lake (Martinez and Anderson, 2013) suggesting methane is emitted in these aerated waters. A unique aspect to the work presented here is the high-frequency sampling over the warm season in Rock Creek which measured consistently high concentrations (400–1,000 nM) in the bottom water (Figure 11) that rival what has been measured in the anoxic bottom waters of the mainstem Chesapeake Bay in mid-summer (Gelesh et al., 2016) and further supports a sustained methane flux to the atmosphere. Thus this relatively shallow (∼3 m) eutrophic estuary, may contribute more methane than previously thought, as was the case for streams and rivers (Stanley et al., 2016), and conforms to our understanding of coastal ecosystems as having an outsized influence on methane fluxes in a global context.

Sedimentary methanogenesis is likely the main source of methane to the water column of Rock Creek because the highest dissolved methane concentrations were measured in the sediments, and methane concentrations in the bottom water were typically higher than the surface water. Biogenic methane is also supported with the sedimentary porewater methane δ¹³C-CH₄ values in the deep sediments being < −70‰ (Whiticar, 1999). Yet, there was also evidence that methane produced in the deep sediments went through some degree of microbial oxidation before reaching the overlying water. The sediment porewater methane profiles showed classic concave-up shapes which are indicative of the anaerobic oxidation of methane (AOM) working in concert with sulfate reduction, as expressed here with the sulfate methane transition (SMT) depths (Jørgensen et al., 2019). AOM is also supported with the porewater methane isotopic composition data. In 2019, the porewater δ¹³C-CH₄ values also increased up the core through the SMT depth. This pattern indicates AOM; as the methane diffuses along the concentration gradient, microbial communities preferentially utilize ¹²C and leave the ¹³C behind (thereby values increase) as methane is oxidized (Whiticar, 1999). In 2019, this pattern is clear regardless of aeration status but there is a shift to more ¹³C depleted values when aerators were OFF. This shift could represent enhanced microbial methane production when aerators were OFF but would need to be validated with other information such as diagenetic modeling (e.g., Martens et al., 1998). The 2019 sedimentary profiles, measured under the low-turbulence diffuse system, support our classic understanding of biogeochemical zonation (Froelich et al., 1979) and suggests that the sediments are diffusion dominated.

The pattern in 2018 was not as clear, possibly due to sediment disturbance with vigorous aeration. The destratification system employed during 2018 involved a high-volume through flow of air that leads to substantial physical disturbance of the water-column and sediments. Such disturbance at the aerators translated into variable porewater methane concentration and the isotope patterns, compared to outside the aeration zone (Figure 5). This implies that a simple, steady-state 1D, diffusion dominated interpretation of these profiles cannot be applied here because the destratification system could have driven substantial advective exchange between sediments and the water-column. It is interesting to note that at the aerators (RC1 and RC2), there was a depletion of ¹³C up the core to the sediment-water interface. One way to inject such a depleted signature is by methanogenesis which could happen in this agitated system by bubbling out deep methane during aeration or methanogenesis in the surface layers using non-competitive substrates (Alperin et al., 1988). More work would be needed to support or refute these possibilities.

Methane from the sediments can either diffuse into the overlying water or bubble out via ebullition (methane oversaturated porewaters forming bubbles). The methane released from both of those processes would have different isotopic signatures. For example, for bubbles to form in the sediments and efflux, the methane concentration must be above saturation which occurs in the sediments deeper than ∼10 cmbsf, depending on the site and year (see Figures 5, 6). At these depths, the δ¹³C-CH₄ values were between −80 and −85‰, so we would expect to see these values in the overlying water column (assuming a small amount of methane from the bubble equilibrates with the water). Bubbles captured in a shallow water column off North Carolina showed no isotopic fractionation when released to the water (Chanton and Martens, 1988). However, if the methane is diffusing into the water from the sediment surface, that methane would carry a δ¹³C-CH₄ value similar to the surface sediment methane signature. For 2018, the δ¹³C-CH₄ values of water column methane were always higher than the surface sediment, regardless of site or aeration status. This suggests that there could be some methane oxidation at the very sediment surface that we are missing in our sediment measurements. This was the same situation when the aerators were OFF in 2019. Perhaps future work could focus on the sediment-water interface as an area of intense methane oxidation, be it either aerobic or anaerobic.

The situation is a little different when the aerators were ON in 2019; the water column methane always had a lower δ¹³C-CH₄ value than the surficial sediments. This suggests that the methane was not simply diffusing in from the sediments, and that there must be another source of methane injecting isotopically depleted carbon; such as, sediment ebullition, methane being advected out of the sediments into the water column from the aerators, or microbial methane production in the water column. We don’t have data to support or refute the advective release of methane out of the sediments other than to say that in 2019, the aeration was much finer and thus most likely less advective than in 2018 so it seems unlikely. Plus, since we were not focused on determining the ebullitive flux, its hard to evaluate. However, in 2018, we had evidence of an ebullitive flux of methane from the sediments. First, the methane in the air above Rock Creek waters had an average δ¹³C-CH₄ value of −50‰ which is slightly depleted in ¹³C from the global, well-mixed value from the northern hemisphere of −47.4‰ (Lan et al., 2021). When this air δ¹³C-CH₄ value was compared to what was measured in the deep sediments, < −70‰, or in the water column, < −60‰, the ¹³C depleted methane in air over Rock Creek waters can be explained with a small amount of biogenic methane from the sediments. Secondly, we measured a direct pulse of biogenic (δ¹³C-CH₄ = −80‰) methane in the air above RC7 at dawn when the aerators were OFF (Figure 7). This was most likely due to ebullition from the sediments directly reaching the atmosphere which was trapped in the air above the water due to the air inversion that occurs at dawn (Crill et al., 1988; Mukhophadhya et al., 2001). While this observation of high concentrations of methane in the morning air is not unexpected, it clearly documents how aqueous environments contribute to the atmospheric methane isotopic signal which has recently been shown to be enhanced in biogenic methane sources (Schaefer et al., 2016; Lan et al., 2021). It also clearly shows that in order to fully capture methane dynamics in this system, future work should quantify the bubble flux from the tributary, much like was done in the temperature lake in California (Martinez and Anderson, 2013).

While we might have captured methane bubbling out of Rock Creek in 2018, that process is very heterogeneous and not continuous, and may not fully explain the water column δ¹³C-CH₄ values from 2019 ON status. We thus explore the possibility of water column methane production under aeration. There is growing evidence for an oxic production pathway in surface waters (see review in Bižić et al., 2020) that might be important. Most notably, it was concluded that 90% of methane in surface waters of a temperate lake was formed in the oxic surface waters, and not the sediments (Donis et al., 2017). However, a recent reevaluation of this work concluded that a sedimentary methane flux from the sediments flanking the lake in the shallow waters could explain the oversaturated surface waters (Peeters et al., 2019). In well stratified lakes, methane production in surface waters is shown to scale with sediment area and mixed layer volume (Günthel et al., 2019). Lakes are hydrodynamically very different from estuaries, making study comparisons to Rock Creek tenuous. Furthermore, the hydrodynamic conditions in destratification systems, such as in 2018, can create large cells of overturning water that may impact a much larger area of the tidal system (Gibbs and Howard-Williams, 2018). The lack of studies of sub-estuaries alone highlights the need to quantify the sources of methane from systems such as Rock Creek.

The box model, which balances sedimentary and air-water methane fluxes, shows that for most of the time, regardless of if aerators were ON or OFF, there is an additional methane source in the water column in 2019 (Figure 11). The observation that outside the aeration zone, RC7 and RC9, there is also an additional methane source could indicate the transport of methane from the Patapsco River, lateral inputs of water from across the creek as destratification cells pull in sources from a cross-section of the creek or simply upstream methane sources. Previous work in rivers also documented higher methane flux in upstream surface waters (de Angelis and Scranton, 1993; Abril and Iversen, 2002; Middelburg et al., 2002); as well as higher methane in smaller width creeks (Borges and Abril, 2011).

The further observation that there is a methane source regardless of the status of the aerators (either ON or OFF) suggests there is a bigger ecosystem response than just the influx of oxygen from the aerators. One possibility could be the presence of algal blooms which produce methane as a byproduct (Bižić et al., 2020). Such blooms have been shown to occur before when aerators are turned OFF (Harris et al., 2015), and we also documented one right after the aerators were turned OFF in 2018 (Figure 10). Increased surface water oxygen levels were observed with an increase in particulate organic nitrogen and carbon, and iron (data not shown). The OsmoSamplers also captured a pulse of methane right after this oxygen pulse. If these blooms are a typical phenomenon right after turning the aerators OFF, it could explain why the methane fluxes were higher 3 days after turning aerators OFF in 2021 (Figure 4D) than when the aerators were ON. Future work could also directly quantify rates of methane production in aerobic waters.

While it seems clear that the aerators are enhancing a methane flux to the air above the creek, that the sediments are the main source of methane to the waters and there possibly is methane production in these aerated waters, we also considered if there is any microbially mediated oxidation in the surface waters that might offset this atmospheric flux. The idea here is that the waters are being oxygenated due to the aeration, as the oxygen data suggested happened during the experiment (Figure 3), and this oxygen is allowing for aerobic methane oxidation to happen. We have already speculated about a small potential oxidation process altering the δ¹³C-CH₄ from our surficial sedimentary methane values to the bottom water. Here, we focus solely on the water column. Since we didn’t directly measure microbial rates, we need to rely on geochemical data. We did this in four ways. First, we interrogated the stable carbon isotopes of methane measured in the water column. The δ¹³C-CH₄ values give a bulk measure of what has happened to that methane since it was formed; the bulk methane signature would be more enriched in ¹³C if it had been oxidized. In order to do this, we considered that the methane in the bottom water is the source for the methane in the surface water. Figures 5, 6 water column isotope profiles clearly show that in all cases, the surface water is always slightly ¹³C depleted, if at all different, contrary to what would be found if aerobic methane oxidation was occurring. Secondly, we looked for a correlation between oxygen concentrations and methane concentrations. For this correlation, we surmised that if oxygen was controlling methane levels, we would see a linear relationship between the two. There was no correlation (Supplementary Figure S8A), which suggests oxygen is not the limiting factor for aerobic methane oxidation. Thirdly, we considered that maybe aerobic methane oxidation was limited by dissolved inorganic nitrogen (DIN), as has been proposed in the literature (Sansone and Martens, 1978). We surmise that if there is a negative relationship between DIN and water column methane concentrations, there could be evidence for aerobic methane oxidation. In other words, if methanotrophs were stimulated by DIN in an otherwise oxic water column, we would find low methane concentrations. This pattern did not emerge (Supplementary Figure S8B). And finally, in the box model approach described in Section 4.3, for most of the time in 2018 and 2019, there was a source of methane and not a sink. The exception was in 2019 when the aerators were ON. Taking all of this geochemical evidence into consideration, we conclude that we found little to no geochemical evidence for water column aerobic methane oxidation. Future work could focus on looking for microbial signature of oxidation in the water column.

We have discussed the differences in atmospheric methane flux between the 2018 and 2019 aeration experiments to most likely be a factor of the bubble size during aeration. The large bubble aerators resulted in a larger flux of methane to the atmosphere than the small bubble aerators. However there were also differences between the 2 years in how long the aerators were turned OFF. Despite variability in the concentrations, we found almost no effect of aeration when the aerators were turned OFF for a few days (2016, 2021), but there was modest methane accumulation as oxygen was depleted for 3–7 days in 2018 and significant methane increases after 13 days in 2019. Thus, it appears that methane accumulation requires more than a day to emerge after deoxygenation, despite prior studies in Arctic lake systems (McIntosh Marcek et al., 2021) and in the Chesapeake Bay mainstem (Gelesh et al., 2016) that appear to indicate that methane is immediately released from sediments following deoxygenation. However, the temporal resolution of prior measures likely cannot capture day to day dynamics, and there are few, if any real-time measures of methane and oxygen in coastal systems to confirm the time-scale of methane responses to deoxygenation.

A second source of complication in our interpretations is that the experimental study is embedded in a system with background environmental variability. Thus, our “ON-OFF” study design using whole ecosystem manipulations does not offer a true controlled experimental system. For example, unlike previous work in Rock Creek (Harris et al., 2015), the waters did not go anoxic within 1 day in 2018 and 2019 and only at RC2 in 2018 did anoxia ever emerge. In 2018, there was also an influx of oxygen-rich water that appeared in the aeration zone within days after the aerators were turned OFF and when oxygen should have been depleted (seen in time series data on Figure 10). Although we did not measure methane during this event, the rapid increase in methane concentrations immediately after this pulse of tidally driven, oxygen-rich water could be the result of a rapid response to deoxygenation (that was not interrupted by the event) or an influx of methane rich water from upstream after the event. Our measurements in 2021 did indicate higher methane concentrations upstream of the aeration zone (∼1,000 nM; Figure 4D), but these are substantially smaller than the concentrations measured after the event in 2018 (1,500 to 2,500 nM; Figure 10), suggesting that an upstream source is unlikely. Other potential inputs of methane are currently unknown, such as groundwater, and there are no substantial tidal wetlands that could serve as a methane source. Although groundwater remains unlikely since chloride concentrations slightly increased with depth in sediment porewaters, our understanding of the drivers for the patterns remains elusive and long term observatories are needed.

Ecosystem responses to aeration beyond deoxygenation also likely feedback to influence methane dynamics in the estuary. Harris et al. (2015) reported a substantial algal bloom in Rock Creek within a week of the aerators being turned off, generating high surface water oxygen and organic matter concentrations. Given that Rock Creek has substantial light attenuation (Secchi Depth < 0.5), the vertical mixing induced by aerators likely keeps phytoplankton mixed below the photic layer, and when aerators are turned OFF, water-column stabilization allows for phytoplankton blooms. This is consistent with recent work in Chesapeake Bay that suggest that high chlorophyll-a packaging under low-light conditions combined with high nutrient concentrations allows for rapid algal growth when light becomes available (Buchanan, 2020). Methane concentrations reached the highest levels we measured at RC2 in 2018 (∼2,700 nM), 3-days after an increase in particulate organic carbon of 700 μM after the aerators were turned OFF (data not shown). The consumption of this organic matter that likely followed may have generated substantial new methane, especially considering that oxygen concentrations were consistently below 32 μM (1 mg L⁻¹) for much of the following week (Figure 10). Given that the “destratification” approach used in 2018 was designed to physically mix and overturn the waters, such a phytoplankton response is likely. The “diffuse” approach used in 2019 involved much less physical disruption of the water-column (water-column salinity profiles were similar during ON and OFF), which may have prevented a phytoplankton response and also allowed for stable conditions that allowed the microbial communities to organize along expected redox conditions. This finding also points to an important impact of aerator design on our results, with differences in impact dependent on the mechanism of aeration and questions remain regarding how best to quantify mixing, oxygenation, and even water fluxes dependent on whether engineering impacted circulation (destratification design through large bubbles) versus diffusion of oxygen (small bubble aeration design).