High-frequency continuous measurements reveal strong diel and seasonal cycling of pCO₂ and CO₂ flux in a mesohaline reach of the Chesapeake Bay

We estimated hourly air–water gas transfer velocities (k₆₀₀) for carbon dioxide in the Rhode River, a mesohaline sub-estuary of the Chesapeake Bay. Gas transfer velocities were calculated from estuary-specific parameterizations developed explicitly for shallow microtidal estuaries in the mid-Atlantic region of the United States, using standardized wind speed measurements. Combining the gas transfer velocity with continuous measurements of pCO₂ in the water and in the overlying atmosphere, we determined the direction and magnitude of CO₂ flux at hourly intervals across a 3-year record (1 July 2018 to 1 July 2021). Continuous year-round measurements enabled us to document strong seasonal cycling, whereby the Rhode River is primarily autotrophic during cold-water months (December–May) and largely net heterotrophic in warm-water months (June–November). Although there is inter-annual variability in CO₂ flux in the Rhode River, the annual mean condition is near carbon neutral. Measurement at high temporal resolution across multiple years revealed that CO₂ flux and apparent trophic status can reverse during a single 24 h period. pCO₂ and CO₂ flux are mediated by temperature effects on biological activity and are inverse to temperature-dependent physical solubility of CO₂ in water. Biological/biogeochemical carbon fixation and mineralization are rapid and extensive, so sufficient sampling frequency is crucial to capture unbiased extremes and central tendencies of these estuarine ecosystems.

1 Introduction

Understanding the air–sea exchange of gases and establishing methodologies for accurate measurements has been a decades-long focus of atmospheric scientists, oceanographers, and biogeochemists seeking to understand interactions between oceans and the atmosphere and how these interactions contribute to the global carbon cycle (Broecker et al., 1979; Wanninkhof, 1992; Wanninkhof et al., 2013). Coastal oceans and estuaries are ecosystems of interest for understanding the complex nature and contribution of the land–sea interface to lateral mass transport of carbon (Abril and Borges, 2005; Cai and Wang, 1998; Frankignoulle et al., 1998; Song et al., 2023) but also with respect to the role these ecosystems play as both atmospheric CO₂ sources and sinks (Abril and Borges, 2005; Chen et al., 2020; Dai et al., 2022; Jiang et al., 2008). The exchange of carbon dioxide, methane, and other greenhouse gases (GHGs) between Earth's atmosphere and inland waters, estuaries, and coastal oceans are well-documented but not fully quantified (Abril and Borges, 2005; Cai, 2011; Laruelle et al., 2017; Raymond and Cole, 2001; Raymond et al., 2013; Van Dam et al., 2019). CO₂ evasion from estuaries alone has been estimated at 15 %–17 % of the total CO₂ input from oceans to the atmosphere (Chen et al., 2020; Laruelle et al., 2017), indicating the regional and global significance of estuaries (Bauer et al., 2013; Frankignoulle et al., 1998; Jiang et al., 2008). Yet, there is still great uncertainty surrounding the true net contributions of coastal oceans, estuaries, and inland water bodies to the atmospheric loading of GHGs (Borges, 2005; Chen et al., 2020; Herrmann et al., 2020; Joesoef et al., 2015; Laruelle et al., 2017; Raymond et al., 2013; Van Dam et al., 2019).

To better understand the effects of estuaries on atmospheric GHG exchange and accumulation, it is imperative that we understand their capacity and function as carbon sources and sinks and ultimately how estuaries factor into the planet's overall global carbon budget (Herrmann et al., 2020; Laruelle et al., 2017; Van Dam et al., 2019). Many attempts to characterize CO₂ flux in estuaries and nearshore oceans (Chen et al., 2013; Herrmann et al., 2020; Rosentreter et al., 2021) have relied on direct measurements using floating domes, tracer gases or, more recently, eddy covariance methods (Laruelle et al., 2017; Van Dam et al., 2019). Because flux measurements are time intensive, they tend to be temporally and spatially limited (Herrmann et al., 2020; Klaus and Vachon, 2020). Using direct flux measurements to derive accurate gas transfer velocity constants (k, the velocity of gas crossing the air–water boundary) enables models to be parameterized to estimate k and compute gas flux. Thus, correlative models that incorporate simultaneous environmental measurements such as wind and/or water velocity, factors that affect turbulence at the air–water interface and promote gas exchange, have aided in the widespread accumulation of gas flux estimates (Raymond and Cole, 2001; Van Dam et al., 2019; Wanninkhof, 2014). Gas transfer velocity constant models vary according to the habitat/system being observed and chemical, physical, and biological factors present in each (e.g., lakes, rivers/streams, estuaries, and oceans; Herrmann et al., 2020; Ho et al., 2016; Raymond and Cole, 2001; Van Dam et al., 2019; Wanninkhof, 1992). To reduce uncertainty in computed gas fluxes, it is critical that the appropriate k models are matched to a targeted ecosystem.

Coastal oceans and estuaries are exceptionally complex, frequently characterized by their relative shallowness and by how their freshwater inputs (riverine, surface, and groundwater) mix with saltwater (Chen et al., 2020). High nutrient and pollutant loading due to urbanization and eutrophication by humans also have important effects on estuaries and coastal oceans (Freeman et al., 2019). High spatial and temporal variability are hallmarks of estuaries.

Here we present a 3-year data set that combines high-frequency (1 min interval) measurements of dissolved and atmospheric CO₂ with co-located and continuous measurements of salinity, water temperature, tidal cycling, and wind velocity, recorded at the Smithsonian Environmental Research Center (SERC) dock in the Rhode River, Maryland. To estimate hourly, daily, seasonal, and annual CO₂ flux rates, we applied a CO₂ gas velocity constant model developed by Van Dam et al. (2019) for the New River, North Carolina. This model is expressly designed for application to shallow well-mixed microtidal estuaries located on the mid-Atlantic coast of the United States.

2 Methods

2.1 Study location

The Rhode River is a tributary and sub-estuary of the Chesapeake Bay, a drowned river valley coastal-plain estuary (Fig. 1). The Rhode River has been studied extensively by SERC staff and colleagues for over 4 decades: nutrient chemistry (Jordan and Correll, 1991; Jordan et al., 1991), phytoplankton ecology (Gallegos et al., 2010), color dissolved organic matter distribution (Tzortziou et al., 2008, 2011), and more recently, modeling of dissolved organic carbon (DOC) input from freshwater and tidal marsh sources have all been researched (Clark et al., 2020). Located on the bay's northwestern shore (38°52^′ N, 76°32^′ W), the Rhode River is bounded at its head by Muddy Creek and at its mouth by the mainstem of the Chesapeake Bay. The Rhode River is a shallow (mean depth = 2 m, max depth = 4.1 m) mesohaline (0 to 18 ppt) well-mixed eutrophic tributary with a length of approximately 5 km; its surface area is approximately 5 km², with a shoreline perimeter of 39 km (Breitburg et al., 2008; Clark et al., 2018). A 0.21 km² tidal marsh (Kirkpatrick Marsh) fringes the estuary at the mouth of Muddy Creek (Fig. 1). Tides are semi-diurnal with a mean amplitude of approximately 30 cm, but water height can be strongly affected by wind and weather events. Muddy Creek is the main freshwater source of the Rhode River and has a maximum flow rate of 10.42 m³ s⁻¹ and mean flow rate 0.18 m³ s⁻¹ (mean flow = 15 552 m³ d⁻¹; Clark et al., 2020, 2018; Jordan et al., 1986). The mean daily volume of freshwater inflow from Muddy Creek is approximately 0.5 % of the mean daily tidal exchange volume based on the Rhode River's area and mean tidal amplitude. In the absence of measurements of the pH or pCO₂ of the freshwater entering the Rhode River from Muddy Creek or other lesser freshwater inputs to the estuary, we are unable to report these pCO₂ or pH values. However, given the exceedingly small overall volume of freshwater input to the Rhode River from its surrounding watershed, it is not considered a river-dominated estuary so is not expected to be substantially influenced by the chemical characteristics of this input. This is not to say that there is no freshwater influence, only that such influences are likely quite local when mixing with far larger volumes of water from the Chesapeake Bay and therefore beyond the resolution of this study.

Although the Rhode River is a model ecosystem that has been studied intensively for several decades across many dimensions (Clark et al., 2018; Correll et al., 1992; Gallegos et al., 1992; Jordan et al., 1991; Rose et al., 2019), no work to date has expressly characterized the nature and dynamics of CO₂ flux between the river and the atmosphere.

Figure 1 Location of study site on the Rhode River, Edgewater, MD (a), as situated in the Chesapeake Bay (b). All pCO₂ and related water quality values reported were measured from the SERC dock, which extends approximately 75 m from the shore into the Rhode River. Red circles indicate locations of the dock and a tidal creek that drains the Kirkpatrick salt marsh (marsh area is equal to 0.21 km², 1 km up estuary from the dock).

2.2 In situ measurements, calculated parameters, and quantities

Continuous, automated environmental measurements were made in and above the Rhode River during a 3-year period between 1 July 2018 and 1 July 2021. The purpose of these measurements was to document fluctuations in aqueous pCO₂ on a fine timescale, from which CO₂ flux between the water and atmosphere could be calculated.

2.2.1 Aqueous CO₂ (pCO_2 water)

To measure the CO₂ gradient ($\mathrm{\Delta }C=p{\mathrm{CO}}_{\mathrm{2}\mathrm{water}}-p{\mathrm{CO}}_{\mathrm{2}\mathrm{air}}$) across the Rhode River surface waters and its overlying atmosphere, measurements of pCO₂ were made with a nondispersive infrared (NDIR) detector. In the case of dissolved gas measurements, water was equilibrated continuously with a spherical falling film equilibrator (Miller et al., 2019). Water from 1 m below the water's surface was pumped and dispersed continuously over a 25.4 cm diameter sphere. The falling film created on the sphere generates a gas exchange surface, which forces CO₂ in the equilibrator headspace into equilibrium with the water's CO₂ content, i.e., mole fraction = xCO₂ (µmol mol⁻¹). Water exits the equilibrator via an airtight drain that prevents headspace contamination from surrounding atmospheric air. Headspace gas circulates continuously in a closed loop through the equilibrator, water trap, and gas dehumidifier; past the NDIR; and back into the equilibrator. Experimental observations concluded that spherical falling film equilibrators achieve 99 % equilibration of CO₂ within 10–15 min, depending on whether step changes are from low to high or high to low; details of the operation and performance of the falling film equilibrator are described in Miller et al. (2019). Measurements were made at 1 min intervals at a pressure equal to the ambient barometric pressure.

Measured raw CO₂ mole fractions (µmol mol⁻¹) were converted to partial pressures (µatm) using Eq. (1). Minute-over-minute values were rounded down to the nearest hour and averaged to provide hourly means. The mole fractions were then evaluated with corresponding water temperature and salinity measurements following the methodology of Zeebe and Wolf-Gladrow (2001), where the saturation vapor pressure of water is calculated according to Weiss and Price (1980) to determine pCO_2 water.

$$\begin{array}{}\text{(1)}& p{\mathrm{CO}}_{\mathrm{2}\mathrm{water}}=x{\mathrm{CO}}_{\mathrm{2}}\cdot (p-p{\mathrm{H}}_{\mathrm{2}}\mathrm{O}),\end{array}$$

where pCO₂ is the partial pressure of CO₂ of water (µatm); xCO₂ is the mole fraction of CO₂ in water (µmol mol⁻¹); p is total pressure, which is equal to 1 atm; and pH₂O is the saturation vapor pressure of water (µatm).

2.2.2 Atmospheric CO₂

Every 6 h, the sample gas stream was automatically diverted with programmed solenoid control valves from the equilibrator to an atmospheric port located approximately 5 m above the pier deck. During atmospheric sampling, 15 1 min interval measurements were made. To account for inaccuracies during the transition period from equilibrator to atmospheric sampling, the final eight measurements were averaged and the first seven were discarded. Similarly, the first 30 measurements following switchover from atmospheric port to equilibrator were discarded, to ensure that measurements were fully equilibrated with water. For these atmospheric measurements, the contribution of the vapor pressure of water to the total atmospheric pressure of the open-air environment was considered negligible (i.e., pH₂O=0 and p=1), such that pCO_2 atm=xCO_2 atm. As such, any potential differences are expected to fall well within the measurement accuracy of the instrument (see below).

One advantage to using a shared NDIR sensor for aquatic and atmospheric samples is that any minor effects of instrument drift will be reflected in both data streams, as opposed to two sensors that drift independently of one another. Likewise, significant and sustained deviation from typical local atmospheric variability will be captured during atmospheric sampling and can signal the need for recalibration and assist with the quality assurance or quality control of corresponding data from both streams. A disadvantage of using a common sensor for both dissolved and atmospheric CO₂ measurements is that it results in a mismatch in sampling frequency of the two. With this limitation in mind, we chose a higher sampling frequency for aquatic measurements to better describe the inherently higher variability in dissolved CO₂ in water vs. that in the atmosphere (Fig. 2).

Figure 2 Hourly pCO_2 water (blue) and pCO_2 air (goldenrod) values from 1 July 2018 to 1 July 2021 (a). The air–water CO₂ gradient, $\mathrm{\Delta }C=p{\mathrm{CO}}_{\mathrm{2}\mathrm{water}}-p{\mathrm{CO}}_{\mathrm{2}\mathrm{air}}$, describes the directionality of gas diffusion. Negative ΔC values (pCO_2 water values falling below the goldenrod demarcation) represent gas movement from air to water and vice versa. Panels (b) and (c) depict mean pCO_2 water values (95 % CI) for each hour of the day for warm and cold seasons, with the dashed lines equal to the mean 3-year value of pCO_2 air.

Given the 3-year time series and strong diel cycling of pCO_2 water (and dissolved oxygen, DO; see Figs. S1 and S2 in the Supplement) in the Rhode River, we chose to aggregate aqueous minute-over-minute measurements to mean hourly averages. Owing to the relative lack of short-term variability in local atmospheric CO₂ concentrations (Fig. 2), we used linear interpolation to impute atmospheric CO₂ concentrations during the hours in between actual readings (6 h gaps between atmospheric measurements), which we assumed to be more realistic and reliable than last observation carried forward (LOCF) methods, where the last observation is repeated for all gaps until the next measurement is encountered, a method that has fallen out of favor, especially for environmental time series data (Lachin, 2016). To determine if any inadvertent bias was introduced by the linear interpolation procedure, summary statistics of actual atmospheric readings and actual readings plus imputed CO₂ values were compared statistically. This approach enabled us to take advantage of >25 000 time points throughout the 3-year period of observation, providing hourly resolution. Mean pCO_2 air=437 ± 20.0 µatm (Table 1), variability that falls well within manufacturer's specifications (see Sect. 2.2.4).

Table 1 Descriptive statistics comparing the seasonality of pCO₂, CO₂ flux, and associated parameters in cold-water (December–May) and warm-water seasons (June–November).

2.2.3 CO₂ gradient (ΔC)

ΔC was determined by subtraction, pCO_2 water−pCO_2 air, where positive ΔC values correspond to higher CO₂ concentrations in the water, tending toward movement from water to air (outgassing or evasion, where Rhode River is a CO₂ source), and negative values signal CO₂ transport from air to water (dissolution, where Rhode River is a CO₂ sink). Figure 2 shows pCO_2 water and pCO_2 air plotted on an hourly basis for the 3-year period beginning 1 July 2018 and ending 1 July 2021. Across this period, ΔC was predominantly negative during cold months and predominantly positive during warm months when pCO_2 water tended to reach the highest values of the year, but ΔC sometimes reversed sign due to occasional extreme daytime photosynthetic drawdown of CO₂ (Fig. 2).

2.2.4 Accuracy of CO₂ measurements

Estimated accuracy of the spherical falling film equilibrator and NDIR sensor (Senseair K30, https://senseair.com/, last access: 15 December 2023) combination was experimentally determined in the lab and found to measure water equilibrated with known gas concentrations to within the ±1 % uncertainty limits of the of certified standard gas mixtures used, well within the published accuracy specification of the Senseair K30 (i.e., ±30 ppmv ± 3 % of instrument reading). Experimental analysis by Martin et al. (2017) report even higher accuracy when relative humidity and atmospheric pressure are controlled for. Details on performance of the spherical falling film equilibrator, such as accuracy, precision, and time constants, can be found in Miller et al. (2019). Although Senseair offers automated calibration via long-term comparisons to atmospheric readings, this feature was deactivated. The K30 NDIR was periodically validated using standard zero CO₂ (nitrogen) and standard certified span gases at intervals of 1–2 months during the study period. Although the K30 was never observed to drift beyond its factory specifications, the sensor was occasionally re-calibrated in the lab, and measured values were accepted without adjustment.

CO₂ measurements were downloaded to a database at approximately 2-week intervals during the observation period. Data were graphed and reviewed visually in combination with twice-weekly observations of equilibrator function recorded in an accompanying notebook. Anomalous data were flagged and excluded from data analysis (e.g., flooding or clogging events that interrupted proper equilibration.)

2.3 Co-located water quality and atmospheric measurements

This water quality station at the SERC dock is a long-term node of the Maryland Department of Natural Resources (DNR) “Eyes on the Bay” Chesapeake Bay tidal water monitoring program and has been operated by the SERC since 1986. Water quality and atmospheric data are maintained by the MarineGEO Upper Chesapeake Bay Observatory and can be accessed online (Benson et al., 2023). A YSI EXO2 sonde was positioned 1 m below the water's surface and in proximity (∼2.5 m distance) to the submerged water pump that fed the pCO₂ equilibrator. Sonde measurements were made at 6 min intervals and aggregated to 1 h averages. The published accuracy specifications for the YSI sonde are as follows: temperature – ±0.01 °C (−5 to 35 °C), salinity – ±1 % of reading or 0.1 ppt (0–70 ppt), and dissolved oxygen – ±0.1 mg L⁻¹ or 1 % of reading (0 to 20 mg L⁻¹). Discrete measurements of temperature and salinity were made using a handheld YSI Professional Plus 2030 with Quattro Cable instrument with the following specifications: temperature – ±0.02 °C (−5 to 70 °C), salinity – ±1 % of reading or 0.1 ppt (0–70 ppt), and dissolved oxygen – ±0.2 mg L⁻¹ or 2 % of reading (0 to 20 mg L⁻¹). Equilibrator temperature was measured with a probe (EDS model OW-TEMP-B3-12xA) accurate to ±0.5 °C (−10 to 85 °C). Discrete measurements were routinely compared to the sonde to corroborate measurement agreement. Wind speed measurements were made using a sonic anemometer (Vaisala WXT-520 weather transmitter) mounted 7 m above the mean low-tide height of the water and located directly above the pCO₂ equilibrator.

2.4 Data processing

Data included in this study span 1 July 2018 to 1 July 2021.

2.4.1 Gas-specific solubility

To determine the purely physical effects of temperature and salinity on CO₂ solubility, gas-specific solubility values K₀ ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}\mathrm{\mu }{\mathrm{atm}}^{-\mathrm{1}}$) were calculated across the 3-year observation period using water temperature and salinity measurements in combination with pCO_2 water values, according to Weiss and Price (1980), at 1 h intervals.

2.4.2 Gas transfer velocity estimation (k)

Given the similarities between the Rhode River and New River estuaries (e.g., shallow microtidal estuaries with slow water velocity and strong diel cycles in pCO₂ and DO), we chose to parameterize gas transfer velocity k (cm h⁻¹) standardized to the unitless Schmidt number 600 (k₆₀₀) according to the estuary-specific k parameterization model developed by Van Dam et al. (2019). Van Dam et al. (2019) determined that k correlated with wind speed differently during daytime vs. nighttime hours (linear vs. parabolic relationships). Wind speed data were collected during the 3-year period from a sonic anemometer located on the SERC dock directly above the equilibration system and approximately 7 m above the water's surface at mean low-tide height. For the analysis, wind speeds were standardized for a height of 10 m following a power-law relationship, $U_{\mathrm{10}}=U_{\mathrm{7}}\cdot (\mathrm{10}/\mathrm{7}{)}^{\mathrm{0.15}}$ (Saucier, 2003). Following Van Dam et al., wind speed data were binned to 1.5 m s⁻¹ intervals for day and night readings and raw values replaced by the mean wind speed for each bin. The median binned wind speed over the Rhode River was 2.2 m s⁻¹ regardless of time of day or season. Recorded wind speeds never exceeded 10 m s⁻¹ and were dominated by much lower values (Fig. S1). Unlike the New River Estuary, the Rhode River's wind speed profile does not differ much between day and night, nor across seasons. For this reason, we chose to use the most conservative k₆₀₀ formulation from Van Dam et al. (2019) that combines day and night winds to estimate k₆₀₀.

Wind speed was used to parameterize k₆₀₀ as follows:

$$\begin{array}{}\text{(2)}& k_{\mathrm{600}}=\mathrm{1.5}\cdot U_{\mathrm{10}}+\mathrm{4.2},\end{array}$$

where U₁₀ is the mean of binned wind speed at 10 m above the water's surface (m s⁻¹).

2.4.3 CO₂ flux

Using continuous, parallel 3-year records (1 July 2018 to 1 July 2021) of dissolved and atmospheric pCO₂, water temperature, salinity, and wind speed (at standard 10 m height; U₁₀), CO₂ flux was derived according to the equation

$$\begin{array}{}\text{(3)}& {\mathrm{CO}}_{\mathrm{2}}\text{flux}=k_{\mathrm{600}}\cdot K_{\mathrm{0}}\cdot \mathrm{\Delta }C\cdot (\mathrm{600}/\mathit{\text{Sc}}{)}^{-\mathrm{0.5}},\end{array}$$

where CO₂ flux is the rate and direction of the CO₂ mass moving between water and gas phases ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$), k₆₀₀ is the gas transfer velocity (cm h⁻¹) normalized to a common Schmidt number (Sc=600), K₀ is gas-specific solubility for CO₂ ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}\mathrm{\mu }{\mathrm{atm}}^{-\mathrm{1}}$), ΔC is the air–water concentration gradient (µatm), and Sc is the Schmidt number.

Note: CO₂ flux calculations require conversion from traditional k₆₀₀ units (cm h⁻¹) to m h⁻¹ and from ΔC units (µatm) to atm prior to calculation.

2.4.4 Day–night designation

To differentiate daytime from nighttime hours, we used the position of the measurements (latitude) in the Rhode River combined with the local date and time. This approach enabled us to uniformly designate various environmental measurements as happening during the day or night (R package LakeMetabolizer; Winslow et al., 2016).

2.4.5 Seasonality

We chose to break the year into two 6-month periods based on seasonal water temperature shifts, designating June–November as “warm-water months” when water temperatures averaged 23.2 ± 6.90 °C (mean ± 1 SD) and December–May as “cold-water months”, 10.9 ± 5.66 °C (Figs. S1 and S2).

2.4.6 Effect size

Owing to the large number of observations available for comparison in this study, the likelihood of finding statistically significant results is quite high. Whether such statistical results by themselves connote practical and informative differences can be difficult to discern. Effect sizes (omega-squared, ω²) were calculated according to two-factor analyses of variance (ANOVAs) where independent variables were investigated by season (cold-water vs. warm-water season) and day–night period and the interaction of season and day–night. The independent variables compared were K₀, CO₂ flux, ΔpCO₂, k₆₀₀, pCO_2 air, pCO_2 water, and wind speed. To account for temporal autocorrelation and lack of independence of observations that are typical of environmental time series data, we corrected for overinflation in the residual mean square used in the effect size calculations by removing the autocorrelation present within residuals, leaving the white-noise component as the unbiased estimate of residual variability (Cochrane–Orcutt procedure, R package orcutt; Spada et al., 2018).

3 Results and discussion

3.1 Daily and seasonal cycling of pCO₂

Hourly averaged measurements of pCO_2 water in the Rhode River across 3 years revealed strong diel and seasonal cycling (Fig. 2). Mean and maximum pCO_2 water were significantly higher in warm-water vs. cold-water months (Table 1). During warm-water months (June–November), daily oscillations of pCO₂ frequently transit from far above to below ambient atmospheric conditions over the course of the day, only to reverse direction (from low to high) during the nighttime hours (Fig. 2). During the summer, pCO_2 water levels sometimes shifted by as much as 4500 µatm in both directions during a single 24 h period (Fig. 2). This pattern is consistent with biologically driven cycling whereby very high early morning pCO_2 water conditions are depleted by net photosynthetic activity (inorganic carbon fixation) over the course of the day, but high pCO_2 water is restored by respiration in the benthos and water column at night (Song et al., 2023). Comparing dissolved oxygen (DO) over the same period, similar harmonic cycling is observed, but maxima and minima of pCO₂ and DO were inversely related (Fig. S1), hallmarks of a production/respiration-driven system (Herrmann et al., 2020; Van Dam et al., 2019).

On the seasonal timescale, pCO₂ was consistently lowest and DO highest during cold-water months of the year (December–May; Fig. S1). Importantly, for both gases the temporal variability (diel cycling; Fig. S2) was most constrained during cold-water months across years, strongly suggesting that carbon fixation exceeds respiration for prolonged periods (weeks to months). In contrast, during warm-water months (June–November), photosynthesis/carbon fixation and respiration are more evenly balanced, compensating for one another over 24 h periods (i.e., respiration > productivity at night and productivity > respiration during daylight hours; Fig. 2).

3.2 Air–water concentration gradient = ΔC (µatm)

When hourly pCO_2 water and pCO_2 air values (composed of 4 hourly measurements and 20 interpolated values per day) were plotted across the 3 years of observation, the diel and seasonal cycles of pCO_2 water are evident. As expected, atmospheric concentrations of CO₂ remained relatively constant compared with aqueous loads. When the mean raw pCO_2 air measurements (mean=435.1, 95 % CI [434.4,435.7]) were compared with raw pCO_2 air measurements plus imputed estimates (mean=435.4, 95 % CI [435.2,435.7]) no statistical difference was observed, indicating that no substantial bias was introduced by linear interpolation of atmospheric measurements.

Although nearshore atmospheric CO₂ concentrations are expected to vary more than those in an isolated well-mixed atmosphere (e.g., at the Mauna Loa Observatory), annual mean values were consistent and within the published uncertainty of the K30 NDIR sensor when compared with global measurements conducted at Mauna Loa (Thoning et al., 2023). Local perturbations (e.g., effects of terrestrial photosynthetic drawdown when wind is absent) were apparent in measurements (Fig. 2), but there were no instances when the measured local atmospheric values were suspiciously high or low for days on end, as compared to expected global mean atmospheric values for the time period (i.e., 408–416 ppmv; https://www.co2.earth/annual-co2, last access: 15 December 2023, Thoning et al., 2023). This lack of sustained anomalous deviation served as additional confirmation that the K30 was functioning properly and had not drifted outside its calibration range. Importantly, given the extreme diel cycling and seasonal variability in the Rhode River's pCO_2 water, the absolute accuracy necessary for determining year-over-year changes in atmospheric or ocean pCO₂ is not a requirement for these CO₂ flux calculations that rely on relative differences between water and atmospheric measurements.

Figure 3 Gas-specific solubility (K₀) for CO₂ based on water temperature and salinity. Units are $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}\mathrm{\mu }{\mathrm{atm}}^{-\mathrm{1}}$ in the Rhode River (1 July 2018 to 1 July 2021).

Hourly air–water concentration gradient values (ΔC = pCO_2 water−pCO_2 air) can be inferred by comparison of the measured concentration values across the 3 years of this study (Fig. 2). During warm months, pCO_2 water routinely shifts from supersaturated to sub-atmospheric and back again over the course of 24 h (e.g., between >2000 and <410 µatm on a single day). These large daily swings in pCO_2 water produced concomitant directional reversals of ΔC (pCO_2 water−pCO_2 air), which result in longer-term averaged gradients (e.g., multi-day or multi-week averages) near zero (Fig. 2). In contrast, most of the time during cold-water months is spent in a state of sub-atmospheric pCO_2 water (under-saturation with respect to the overlying atmosphere), resulting in ΔC values that are negative and that promote movement of CO₂ from the atmosphere into the water.

3.3 Gas-specific solubility (K₀)

To account for the physical effects of temperature and salinity on the solubility of CO₂ in estuarine water, K₀ was calculated using the methods of Weiss and Price (1980). K₀ varied strongly across seasons over the 3-year observation period. The maximum annual range was 0.027 to 0.071 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}\mathrm{\mu }{\mathrm{atm}}^{-\mathrm{1}}$, mean cold-water months were 0.051, and mean warm-water months were 0.035 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}}\mathrm{\mu }{\mathrm{atm}}^{-\mathrm{1}}$, confirming that CO₂ was most soluble during winter and least soluble in summer (Fig. 3). This is inverse to observed dissolved CO₂ values: pCO_2 water was lowest and least variable during winter and highest and most variable during summer (Fig. 2, Table 1), suggesting that solubility in and of itself plays only a minor and non-limiting role in pCO_2 water in the Rhode River. Effect size (ω²) estimates indicated that the greatest proportion of variability in K₀ was associated with season vs. day–night or the interaction of the two (Table 2).

Table 2 Contrast effect sizes based on two-factor ANOVA where independent variables were compared by season (cold-water season, December–May, vs. warm-water season, June–November), day–night period, and the interaction of the two. ω² is a measure of effect size, estimating the proportion of total variance explained by each parameter. Effect sizes were corrected for inherent temporal autocorrelation using the Cochrane–Orcutt procedure (Spada et al., 2018).

3.4 Temperature–biology ratio

To independently parse the magnitude of the physical vs. biological forcing of pCO_2 water, we normalized values to mean water temperature and estimated the Takahashi temperature–biology ratio (Takahashi et al., 2002) to compare the influence of temperature and biological activities on pCO_2 water. Across the 3-year period, we found that the predominant driver of pCO_2 water in the Rhode River was biological activity, accounting for nearly 4 times more forcing than the physical effects of water temperature on CO₂ solubility (Table 3). These patterns demonstrate the outsized role that biological processes play in shaping pCO_2 water in nearshore marine and estuarine ecosystems (Dai et al., 2022; Van Dam et al., 2019).

Table 3 Takahashi temperature–biology ratio (Eq. 5a from Takahashi et al., 2002).

3.5 Gas transfer velocity (k₆₀₀)

Gas transfer velocity is affected by both mass transfer from molecular diffusion driven by ΔC (i.e., CO₂ gradient between water and the atmosphere) and momentum transfer linked to external environmental forces that enhance turbulence at the air–water boundary layer (Ho et al., 2016; Raymond and Cole, 2001; Van Dam et al., 2019). Van Dam et al. (2019) validated the use of wind speed at 10 m above the water's surface (U₁₀) to estimate gas transfer velocities of CO₂ that were standardized to a Schmidt number of 600 (k₆₀₀) by comparing estimated values to k₆₀₀ values derived directly from eddy covariance CO₂ flux measurements. Given the relative uniformity of wind speed over the Rhode River, where median binned U₁₀ wind speed (converted from U₇ measurements) was 2.2 m s⁻¹ regardless of time of day or season and maximum values that rarely exceeded 10 m s⁻¹ (Table 1, Fig. S1), we chose to use the most conservative estuarine-specific parameterization of k₆₀₀ (Van Dam et al., 2019; Eq. 2). The mean overall Rhode River k₆₀₀ value for CO₂ (mean ± SD, 7.86 ± 2.05 cm h⁻¹) was of comparable magnitude to that of the New River Estuary, NC (9.37 ± 9.47 cm h⁻¹). However, wind speed varied far less on the Rhode River than the New River Estuary and day–night explained more variability in wind speed than season. Because wind speed directly influenced the formulation of k₆₀₀ (Eq. 2), the effect size of day–night is similarly greater than the seasonal effect on gas transfer velocity (Table 2). Nevertheless, effect sizes (ω²) indicate that season explained at least 10 times more of the observed variance of pCO_2 water, pCO_2 air, air–water concentration gradient, CO₂ flux, and gas-specific solubility than day–night or their interaction (Table 2). Given the minor freshwater input and microtidal nature of the Rhode River, we do not believe that lateral water velocity and bottom turbulence appreciably affect the gas transfer velocity of CO₂ here, although we did not investigate possible influences explicitly.

Importantly, in coastal marine and estuarine habitats, ΔC can shift as much as several thousand micro-atmospheres during a single day due to diel cycling associated with CO₂ production and depletion (Figs. 2 and S2). The uncertainty surrounding gas transfer velocity parameterization can represent a major source of error in CO₂ flux calculations (Frankignoulle et al., 1998; Upstill-Goddard, 2006; Wanninkhof and McGillis, 1999); however, small errors in k₆₀₀ have far less effect on CO₂ flux calculations in estuaries, which experience pCO₂ swings of several thousand micro-atmospheres during a single day, compared with more stable conditions of the open ocean where interannual ranges of pCO₂ are typically far less (Van Dam et al., 2019).

3.6 CO₂ flux–seasonality and interannual variation

CO₂ flux was determined according to Eq. (3) using hourly ΔC measurements, CO₂ solubility values (K₀) calculated according to temperature and salinity, and the estuary-specific standardized gas transfer velocities (k₆₀₀) of Van Dam et al. (2019). CO₂ flux was plotted across the 3 years of observations at hourly and monthly intervals (Fig. 4a and b). As observed with pCO₂, CO₂ flux in the Rhode River was shown to be strongly seasonal. Given the similarity in wind speed across seasons (Fig. S1), the effect of differential mean ΔC and variation between warm- and cold-water seasons (Fig. 2, Table 1) almost certainly drive the observed seasonal differences in CO₂ flux (Fig. 4). Again, the specific solubility of CO₂ is greatest at low temperatures, yet this is contrary to the observed mean pCO_2 water patterns, pointing toward a biological mechanism for pCO₂, ΔC, and ultimately, CO₂ flux. The effect size of season on CO₂ flux was 2 orders of magnitude greater than either the individual factor of day–night or the interaction of factors, season x day–night (Table 2).

Figure 4 CO₂ flux estimates by year: (a) hourly and (b) monthly average CO₂ flux estimates with 95 % confidence limits. Black dots in panel (b) indicate mean monthly fluxes across years.

Among years, pCO_2 water and CO₂ flux largely repeat themselves, with dissolved CO₂ becoming consistently sub-atmospheric and CO₂ flux going negative (gas exchange from atmosphere to water) between December and May and abruptly transitioning to much higher maximum, yet variable, pCO_2 water values with net positive CO₂ fluxes from June through November (Figs. 2 and 4). Monthly averaged CO₂ fluxes are consistent among years (Fig. 4b), with net positive CO₂ fluxes (heterotrophic conditions) between June and November and negative (autotrophic) fluxes dominating when water temperatures are cold, between December and May. Despite the overall similarities in seasonal CO₂ flux, inter-annual patterns can vary considerably. When hourly CO₂ flux values were averaged for the year, the Rhode River in 2019 was shown to have a net positive flux, but it had a net negative flux in 2020. When scaled for the year, 2019 outgassed CO₂ from the water to the atmosphere at a rate of 2215.08 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (95 % CI = 1816.88, 2613.29). The annual net flux rate in 2020 was negative (i.e., CO₂ moved from the atmosphere into the river) at a rate of −1361.31 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (95 % CI = −1723.60, −999.01).

At shorter timescales, such as comparing the same week of the year among years, we sometimes observed vast differences in the magnitude and direction of CO₂ flux (Fig. S3 in the Supplement), signaling differences in seasonal conditions between years. Transient events can also result in deviations from otherwise typical CO₂ flux conditions. For example, the period from July 2018 to January 2019 deviated from other years as CO₂ flux was more erratic, with intermittent periods of negative and positive CO₂ flux extending later into the winter season than in other years. When water temperatures are compared among years, 2018 was shown to be more inconsistent, with more pronounced temperature shifts and reversals than in 2019 or 2020 (Fig. S1). Salinities remained relatively low for the latter half of 2018 into early 2019, reflecting wetter conditions (Fig. S1). There were also two rapid salinity declines (>4 ppt reductions) in July and October 2018, likely associated with strong precipitation events. These events were both followed by immediate spikes in chlorophyll-a concentration to levels exceeding 200 µg L⁻¹, indicative of phytoplankton bloom conditions. From 2018 to 2021, chlorophyll-a levels of this magnitude and greater were generally confined to cold-water months (December–May; Fig. S1). Erratic water temperature and salinity are also reflected in more variable gas-specific solubility (K₀) for CO₂ in 2018 than in later years (Fig. 3).

Gallegos et al. (1992) documented predictable phytoplankton blooms associated with freshets in the Rhode River, when nutrient-rich freshwater inundates the estuary not from point and non-point sources within the local Rhode River watershed but instead from the enormous watershed that feeds the Susquehanna River, the primary source of freshwater input into the Chesapeake Bay above the Potomac as well as >50 % of the entire bay's freshwater (US Geological Survey, 2023). Unlike river-dominated estuaries, in the Rhode River estuary, volumetric influxes from the Chesapeake Bay end member far exceed freshwater input from the Muddy Creek and secondary tributaries. In the Rhode River, phytoplankton blooms result in the temporary depletion of pCO_2 water, followed by a spike as phytoplankton senesce and organic carbon is decomposed/re-mineralized back into inorganic carbon. Episodic, short-lived occurrences like these demonstrate how immediate small-scale biological forcing can be coupled with and catalyzed by distant large-scale weather and hydrological events. These in turn can influence pCO₂ flux variations within seasons and among years (Figs. 3 and S3 and Chen et al., 2020).

Overall, except for wind speed, the effect sizes for the other six measured or calculated variables were shown to be greatest for season vs. day–night or season vs. the interaction of season x day–night, and in all cases the season effect was greater by at least 1 order of magnitude (Table 2). Seasonality has 10 to 1000 times more explanatory power than other variables investigated, as estimated by ω² (Table 2).

3.7 Shifting net ecosystem production

To better understand how the net ecosystem production (NEP) of the Rhode River shifts throughout the year, where positive NEP indicates that the river is storing carbon (autotrophic state) and negative NEP indicates that it is releasing carbon to the atmosphere (heterotrophic state), we calculated hourly CO₂ flux values, averaged them by day (i.e., 24 h period), and plotted each in relation to the ΔC=0 reference. Each day of the 3-year study was categorized as either net heterotrophic (CO₂ flux from water to atmosphere) or net autotrophic (CO₂ flux from atmosphere to water). Each day was then further identified as purely heterotrophic (all hours of the day were heterotrophic), purely autotrophic, or mixed (some hours were heterotrophic and some were autotrophic, resulting in a net autotrophic or net heterotrophic state for the day; Fig. 5). From July 2018 to July 2021, most 24 h periods were categorized as pure autotrophic ($\mathrm{444}/\mathrm{920}=\mathrm{48}$ %), while 25 % ($\mathrm{229}/\mathrm{990}$) were purely heterotrophic and the remainder of mixed trophic status (17 % net heterotrophic and 10 % net autotrophic; Fig. 5).

Figure 5 Daily mean CO₂ flux estimates (the CO₂ gradient is CO_2 water−CO_2 air). Green dots indicate days when all 24 hourly flux measurements were negative (autotrophic with +NEP); blue dots indicate days on which all 24 hourly flux measurements were positive (heterotrophic with −NEP). Orange dots indicate that hourly fluxes were both negative and positive, and the position of the orange dot below or above the zero line indicates whether the day was net autotrophic or net heterotrophic. The insert describes the proportion of days in each category, indicating that during 58 % (0.48+0.10) of days across 3 years of observation, the Rhode River was a CO₂ sink.

Altogether, the Rhode River was net autotrophic for 535 of 920 d (58 %) and net heterotrophic for 42 % of the time (385 d) across 3 years. When CO₂ flux is integrated over all 3 years, the Rhode River is shown to have near-neutral NEP (Fig. 6). The effect size of season is 2 orders of magnitude greater than either that of day–night or of the season x day–night interaction (Table 2). Mean CO₂ flux values highlight the obvious correlation between season and NEP; error bars (±1 SD) reveal the importance of diel cycling where the magnitude and directionality of day–night flux variability is approximately equal to the overall variability accrued across all 3 years (Fig. 6). Although CO₂ flux is less variable and more autotrophic during cold-water months than warm-water months in the Rhode River, the range of possible values that occur across night and day, regardless of season, must be taken into consideration to minimize incidental sampling bias (Figs. 2 and 6).

Figure 6 Mean CO₂ flux ± 1 SD ($\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) plotted by day–night cycling, cold-water/warm-water season, season vs. day–night interaction, and overall CO₂ flux across 3 years of observation.

A multi-year investigation of CO₂ flux in the mainstem of the Chesapeake Bay by Chen et al. (2020) combined several bay-wide cruises that were distributed across seasons to collect discrete and underway pCO₂ data for CO₂ flux calculations. They concluded that the low-salinity upper bay, which receives large volumes of freshwater from the Susquehanna River, was net heterotrophic; the mesohaline middle bay was net autotrophic, and the polyhaline lower bay was near carbon neutral. Chen et al. (2020) characterized Chesapeake Bay on the whole as a weak source of CO₂ to the atmosphere (net flux = 0.73 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) but suggested that during wet years, it may function as weak sink of CO₂. Herrmann et al. (2020) also concluded that the Chesapeake Bay was a weak source of CO₂ to the atmosphere based on calculated pCO₂ values from long-term pH and alkalinity measurements (net flux = 1.2 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Brodeur et al. (2019) examined dissolved inorganic carbon (DIC) and total alkalinity along the mainstem of the Chesapeake Bay across the year in 2016 and concluded that DIC increases from north to south and from surface waters to depth but that seasonal riverine input and biological cycling were significantly important, concluding that the bay as a whole was a sink for CO₂.

When our annual mean pCO₂ values were compared with the Chen et al. (2020) survey, the Rhode River was shown to be higher on average and more variable than the mesohaline mainstem of the bay (591 ± 652 vs. 416 ± 167 µatm), including a substantially greater measured range (min=15, max=5182 µatm vs. 103 and 1033 µatm). These results suggest that water in the shallow and well-mixed Rhode River, and DIC in particular, undergo more acute biological transformation than in the mesohaline mainstem of Chesapeake Bay. Chen et al. (2020) point to a variety of factors that affect pCO₂ and CO₂ flux in the mainstem bay, including temperature, depth, stratification, and freshwater input volume, some of which may attenuate biological forcing. Interannual variability was demonstrated in both the Rhode River (some years were net autotrophic and others heterotrophic; Figs. 4 and 5) and in the mesohaline mainstem of the bay; however, we attribute interannual variability in pCO₂ and CO₂ flux primarily to variation in water temperature that in turn drives biological activity. We conclude that seasonal variations in the Rhode River (and likely similar rivers in the mesohaline portion of the Chesapeake) are significant and predictable, closely associated with water temperature, and that temperature mediates NEP biologically rather than by changes to the solubility of CO₂.

In the Rhode River, we find that CO₂ flux reverses itself daily for part of the year (June–November) yielding some days that are characterized as a net sink (net autotrophic) and others that are a net source (net heterotrophic). From December to May, diel cycling is minimal and the river is almost exclusively a net sink, autotrophic both day and night. Finally, although CO₂ flux is pronounced but variable across seasons, the net CO₂ flux of the Rhode River on an annual basis is near neutral.

3.8 Lateral transport

Tidal cycling has been shown to liberate and laterally transport DOC from brackish marshes to adjacent estuaries (Cai, 2011; Herrmann et al., 2015) and therefore is of great importance to carbon cycling and budgets of wetlands and estuaries (Najjar et al., 2020). DOC outwelling from the Kirkpatrick Marsh (hereafter KPM), a 0.21 km² tidal marsh located approximately 1 km up estuary from our primary study site at the SERC Dock (Fig. 1), into the Rhode River has been measured and modeled extensively in recent years (Clark et al., 2020; Menendez et al., 2022; Tzortziou et al., 2011, 2008). These studies indicate that the KPM is responsible for a large portion of overall DOC input to the Rhode River, as well as significant export from the river to the mainstem of Chesapeake Bay. Model generation and validation by Clark et al. (2020) indicate that up to 13.1 % of the total DOC input to the Rhode River originates in the KPM. Another important source (53 % of total) is DOC derived from phytoplankton, which is therefore labile and readily biodegraded and remineralized into DIC. Furthermore, large quantities of other semi-labile forms of DOC are exported from the KPM, which are themselves subject to photochemical processes and biodegradation and remineralization (Clark et al., 2020). Importantly, each of these DOC streams provides a potential source of DIC, including pCO₂, to the Rhode River.

Dissolved inorganic carbon generated in brackish tidal wetlands is also outwelled directly into estuaries (e.g., Cai et al., 2000; Chu et al., 2018; Song et al., 2023). Recent work by Song et al. (2023) demonstrates that pCO₂ in a salt marsh tidal creek in Waquoit Bay, MA, was regulated by both tide height (inversely) and the day–night cycle, with nighttime low tides resulting in the highest pCO₂ values, signaling a strong local effect from respiration and photosynthesis in combination with tidal outwelling.

In the Rhode River watershed, pCO₂ was measured continuously in the single tidal creek that drains the KPM using the same methods as at our primary study location. We observed that the KPM tidal creek pCO₂ follows the tidal cycle exclusively, yet outside the mouth of the tidal creek in the estuary proper, day–night cycling overwhelms this marsh tidal signal. Simultaneous pCO₂ measurements from the SERC dock follow a strict day–night cycle (Fig. S4 in the Supplement). However, while peak levels of dissolved CO₂ in the Kirkpatrick Marsh creek occur at low tide and can reach values nearly 20 times greater than highs at the SERC dock (Fig. S4), there is no obvious evidence of this tidal DIC input at the dock site. Remineralization of DOC exported from the KPM, as well as DOC originating in other locations within the watershed, are important sources of DIC in the river; however, given the relative volumes of these sources to that of the much larger estuary as well as the physical distance (∼1 km) from SERC dock, these input signals should be expected to lag and be damped inside the estuary and not tightly coupled with tidal cycles. Instead, pCO₂ exported from the KPM is expected to undergo significant dilution effects, be partially off-gassed to the atmosphere, and be metabolized via photosynthesis, reducing its influence on downstream sites. These findings suggest that despite periodic extreme pCO₂ in KPM tidal creek (>30 000 ppmv), the overall mass of CO₂ export is not sufficient to have measurable effects on the deeper, well-mixed portions of the Rhode River.

Thus, although land–sea interfaces and the outwelling of DOC and DIC are important in estuaries and coastal ecosystems, the relative sizes of wetlands and adjacent water bodies and the overall volume of water moving between the two are also important factors. In eutrophic estuaries like the Rhode River, biological forcing can rapidly assimilate DIC and degrade and mineralize labile forms of DOC, as evidenced by extensive diel cycling in these systems (e.g., Brodeur et al., 2019; Song et al., 2023, and the present study.) The much larger and more complex Chesapeake Bay generally follows seasonal changes in pCO₂ and CO₂ flux, but these appear to be most predictable in the upper oligohaline portion and in the polyhaline region of the bay near the mouth, where freshwater and oceanic end-member effects are most pronounced (Brodeur et al., 2019; Chen et al., 2020). The central mesohaline part of Chesapeake Bay comprises numerous discrete and unique watersheds and sub-estuaries/rivers, each of which exchanges water with the bay. Elucidating spatial and temporal patterns of pCO₂ and CO₂ flux are vital for understanding each one's role as an atmospheric source or sink but could also provide better insight into how each may be influenced by global increases in atmospheric CO₂ (i.e., acidification and its influences on estuarine metabolism and the local biota, fisheries, and habitats each support.) Collectively, these and other sub-estuaries will have cumulative effects on the overall water quality of Chesapeake Bay, including cycling of DOC and DIC, which in turn affect pCO₂ and CO₂ flux.

4 Conclusion

The notion that estuaries are predominantly heterotrophic systems that invariably outgas more CO₂ to the atmosphere than they absorb has been a long-held view (Abril et al., 2000; Borges et al., 2004; Cai, 2011; Cai et al., 2000; Chen et al., 2013; Frankignoulle et al., 1998; Gattuso et al., 1998). However, more recently investigators have realized that physical and hydrological characteristics, geographical location, size, and biological and biogeochemical activities may individually, or together, influence CO₂ flux in estuaries and therefore contributions to atmospheric chemistry (Brodeur et al., 2019; Caffrey, 2004; Chen et al., 2013, 2020; Herrmann et al., 2020). As indicated in this study and others, the role that biological processes play in estuaries to either fix CO₂ (autotrophy) or liberate CO₂ (heterotrophy) are extensive, complex, and can be quite variable over space and time (Brodeur et al., 2019; Chen et al., 2020; Herrmann et al., 2020; Rosentreter et al., 2021). High-frequency automated measurements revealed strong seasonal contrasts in dissolved CO₂ content and CO₂ flux between water and the atmosphere of the Rhode River, a shallow mesohaline reach of the Chesapeake Bay. Importantly, only through high-frequency multi-year measurements could diel and seasonal cycling be fully discerned. Indeed, inadequate sampling can induce bias (e.g., upscaling from a small number of daytime samples taken during warm-water months can skew apparent patterns; Laruelle et al., 2017; Van Dam et al., 2019). The timing and frequency of measurements are critical and have potential for strong and misleading biases if sampling is insufficient. In contrast, cold-water months coincide with long periods (weeks to months) of continuous sub-atmospheric sink conditions for CO₂. Using these measurements, we estimated the direction and magnitude of CO₂ flux in hourly, daily, and annual terms. In the Rhode River, CO₂ flux reverses itself daily for part of the year (June through November), yielding some days that are characterized as net sink (net autotrophic and NEP>0) and others that are net source (net heterotrophic and NEP<0). From December to May, diel cycling is minimal, and the river is almost exclusively a CO₂ sink with +NEP both day and night. Although CO₂ flux is pronounced but variable across seasons, the net CO₂ flux of the Rhode River on an annual basis is near carbon neutral, although some years are net heterotrophic and others net autotrophic.

High-frequency sampling of pCO₂, although typically confined spatially, is one approach to understanding fundamental aspects of estuarine metabolic states and CO₂ flux that may otherwise go undetected (Song et al., 2023). To address the spatial complexity of estuarine, nearshore, and inland waters, more observation locations are required. As with any environmental or ecological question, careful sampling design is critical to balance efficiency and statistical power.

As the largest and arguably most complex estuary in the United States, the Chesapeake Bay is the subject of extensive ecosystem management efforts and ranks among the most studied and monitored estuaries in the world (Boesch and Goldman 2009). Yet, information on CO₂ and GHG fluxes continues to be limited (Brodeur et al., 2019; Chen et al., 2020; Herrmann et al., 2020). Given the extensive coordinated monitoring programs that either make real-time water quality measurements and/or maintain routine water sampling schedules (e.g., the Maryland DNR Eyes on the Bay program) in this region, existing water quality observation assets and sampling programs could be leveraged to more fully characterize and quantify CO₂ and other GHG dynamics and flux in the bay and elsewhere (see Saba et al., 2019). For example, coordinated deployment of additional automated sampling devices (e.g., robust air–water equilibrators and traditional atmospheric gas sensors) in key locations would enable estimates of CO₂ flux, and if combined with pH, DIC, or total alkalinity measurements, carbonate chemistry calculations as well. Importantly, such installations need not be permanent. Instead, a small group of instruments could be systematically deployed across an existing observation network, co-located with other water quality instruments using a stratified sampling approach to capture spatial variability. For example, a set of shifting 2–4-week deployments during summer and winter months could yield sufficient data to advance our understanding of Chesapeake-Bay-wide CO₂ flux significantly in a single year. Such information would complement underway transects that are vital but that tend to underestimate temporal variability in any given location. In the case of dissolved GHGs, liquid–air equilibration techniques are being used to measure multiple GHGs (Call et al., 2015; Hartmann et al., 2018; Gülzow et al., 2011; Miller et al., 2019; Xiao et al., 2020).

Understanding the GHG dynamics in estuaries is a vital component to generating accurate global budgets (Maher and Eyre, 2012) as well as informing where emerging carbon capture technologies, including nature-based solutions, might be best located (Bradshaw and Dance, 2005; Sun et al., 2021). In the case of estuaries, there have been extensive global losses of seagrasses due to habitat degradation, pollution, and disease (Waycott et al., 2009). In addition to many other ecosystem service benefits, restoration of seagrass and submerged aquatic vegetation has the potential to restore and enhance natural carbon sequestration (i.e., blue carbon; Kennedy et al., 2022; Macreadie et al., 2022; Unsworth et al., 2022). In Virginia, USA, Oreska et al. (2020) demonstrated how the functional benefits of a restored seagrass meadow habitat can be quantified ecologically in terms of its ability to sequester carbon and affect GHG fluxes between the estuary and atmosphere. Uniquely, these investigators then monetized the costs and benefits of habitat restoration and function as CO₂ offset credits as part of a GHG budget and demonstrated how such approaches can be used to incentivize habitat restoration (Oreska et al., 2020).

Increasing the completeness and utility of global GHG budgets as they relate to human activities and ecosystem functions are necessary steps toward combating global climate change. Measurement of GHGs at high spatial and temporal resolution using economical automated measurement solutions can increase our understanding of GHG dynamics at small, ecologically significant scales as well as at the larger ecosystem level of an estuary.