Challenges in Quantifying Air-Water Carbon Dioxide Flux Using Estuarine Water Quality Data: Case Study for Chesapeake Bay.

Estuaries play an uncertain but potentially important role in the global carbon cycle via CO2 outgassing. The uncertainty mainly stems from the paucity of studies that document the full spatial and temporal variability of estuarine surface water partial pressure of carbon dioxide ( pCO2). Here, we explore the potential of utilizing the abundance of pH data from historical water quality monitoring programs to fill the data void via a case study of the mainstem Chesapeake Bay (eastern United States). We calculate pCO2 and the air-water CO2 flux at monthly resolution from 1998 to 2018 from tidal fresh to polyhaline waters, paying special attention to the error estimation. The biggest error is due to the pH measurement error, and errors due to the gas transfer velocity, temporal sampling, the alkalinity mixing model, and the organic alkalinity estimation are 72%, 27%, 15%, and 5%, respectively, of the error due to pH. Seasonal, interannual, and spatial variability in the air-water flux and surface pCO2 is high, and a correlation analysis with oxygen reveals that this variability is driven largely by biological processes. Averaged over 1998-2018, the mainstem bay is a weak net source of CO2 to the atmosphere of 1.2 (1.1, 1.4) mol m-2 yr-1 (best estimate and 95% confidence interval). Our findings suggest that the abundance of historical pH measurements in estuaries around the globe should be mined in order to constrain the large spatial and temporal variability of the CO2 exchange between estuaries and the atmosphere. ©2020. The Authors.

Introduction

Estuaries outgas a significant amount of carbon dioxide (CO2) to the atmosphere, ranging from 0.1 to 0.5 Pg C yr−1, up to one quarter of the uptake of carbon dioxide by the open ocean (Borges & Abril, 2011; Cai, 2011; Chen et al., 2013, 2012; Chen & Borges, 2009; Laruelle et al., 2010). The uncertainty in the outgassing stems from the large spatial and temporal gradients in the partial pressure of CO2 (pCO2) in estuarine surface waters, as high as ~100 μatm km−1 and >1,000 μatm month−1, respectively (Borges & Abril, 2011), combined with limited observations. In the most recent global synthesis of estuarine CO2 outgassing (Chen et al., 2013), data from 165 estuaries were included, but only 29 systems had seasonal coverage and only 16 of these had published studies of their pCO2 dynamics (Àlvarez et al., 1999; Frankignoulle et al., 1998; Guo et al., 2009; Jiang et al., 2008; Koné et al., 2009; Mukhopadhyay et al., 2002; Raymond et al., 2000; Zhai et al., 2007). Fewer still had interannual coverage, with 13 of the 16 published seasonal studies covering a time period of less than 2 yr and the remaining 3 consisting of up to nine cruises over as many as 5 yr. Of the studies published since the synthesis of Chen et al. (2013), we are only aware of three that cover more than a few years (Carstensen et al., 2018; Dinauer & Mucci, 2017; Prasad et al., 2013).

This study is based on the premise that historical water quality data can help fill the temporal and spatial variability void in estuarine surface water pCO2. Support for this premise is found in the multidecadal pCO2 studies of Carstensen et al. (2018) and Prasad et al. (2013), who computed pCO2 using pH and total alkalinity (TA) data from monitoring programs, which increased in number in the 1970s and 1980s in response to a growing public awareness of water pollution. Only recently have monitoring program data sets begun to be used for long‐term studies of the estuarine carbonate system, with a focus on acidification and alkalinization (Baumann & Smith, 2017; Carstensen & Duarte, 2019; Duarte et al., 2013, Hu et al., 2015; Najjar et al., 2020; Waldbusser et al., 2013). However, there are at least three challenges associated with using water quality data for studying pCO2 dynamics. First, water quality monitoring pH and TA data are of relatively low quality when compared to pH and TA measured using modern methods, which translates into a large uncertainty in computed pCO2. For example, for conditions typical of a midlatitude estuary (temperature = 15°C, salinity = 15, pH = 8, and TA = 1.5 mol m−3), the accuracy of pH probes typically used in water quality monitoring programs (0.2) translates to a pCO2 accuracy of 50%. Second, monitoring programs either sample on a roughly monthly time interval throughout a system (e.g., Chesapeake Bay Program, 1996), thereby missing variability due to tides, solar diurnal forcing, and weather events or continuously at a few locations within a system (e.g., Baumann & Smith, 2017), thereby missing the large spatial variability from tidal fresh to polyhaline waters. Third, compared to oceanic waters, the carbonate system in estuaries is poorly constrained because of uncertainties about equilibrium constants and organic alkalinity (Dinauer & Mucci, 2017). To our knowledge, no long‐term studies of surface water pCO2 and the air‐water CO2 flux have addressed these multiple challenges and associated sources of error.

In this study, we present an analysis of the surface water pCO2 and the air‐water CO2 flux for the mainstem Chesapeake Bay using measurements of pH, temperature, and salinity from a water quality monitoring program. The spatial and temporal coverage of the analysis is unprecedented for an estuary, based on monthly or semimonthly surveys of 33 stations throughout the mainstem bay from 1998 to 2018. Special attention is given to error estimation by propagating errors in pH, TA, organic alkalinity, temporal sampling, and the gas transfer velocity. Particularly novel is the use of high‐quality measurements of TA and dissolved inorganic carbon (DIC) to better constrain errors. We document seasonal, interannual, and spatial variability of the computed surface water pCO2 and the air‐water CO2 flux and identify the likely driving factors for this variability through statistical analysis. Dissolved oxygen supersaturation (Δ[O2], the negative of the apparent oxygen utilization) is used as a tracer of net ecosystem production (NEP = photosynthesis minus respiration) and a means for assessing NEP's influence on surface water pCO2.

Study Area, Observations, and Methods

Study Area

The mainstem Chesapeake Bay is a large, partially stratified, coastal plain estuary located in eastern North America (Figure 1). The standard Chesapeake Bay management segmentation scheme (Chesapeake Bay Program, 2004) partitions the bay into eight segments (Table 1 and Figure 1). Segment 1 is considered tidal fresh, Segment 2 oligohaline, Segments 3–5 mesohaline, and Segments 6–8 polyhaline (see salinity ranges in Table 1). The mainstem Chesapeake Bay ranks high in terms of primary productivity among the world's estuaries (Cloern et al., 2014). Kemp et al. (1997) found the bay to be net heterotrophic in Segments 1–3, presumably resulting from a combination of high turbidity that reduces photosynthesis and high riverine organic carbon loads that stimulate respiration. The clearer waters of Segments 4–5 and 6–8 were found to be metabolically neutral and net autotrophic, respectively. Overall, the bay is thought to be net autotrophic (Kemp et al., 1997, C. Shen et al., 2019), which is unusual among the world's estuaries (Borges & Abril, 2011).

Figure 1

Mainstem bay segmentation and locations of the USGS nontidal alkalinity station and tidal pH, temperature, salinity, and dissolved oxygen stations from the Chesapeake bay Program (CBP). Stations that were used for an evaluation of the pH data quality, as described in the supporting information, are circled (CB5.3 and 5.4). Also shown are Shadwick, Friedrichs, et al. (2019) buoy and Shadwick, De Meo, et al. (2019) stations.

Table 1

Characteristics of the Eight Segments of the Mainstem Chesapeake Bay

Segment	Salinity regime	Seasonal salinity range a	Area b (km2)	Mean depth c (m)	Station names
1	Tidal fresh	0–1	152	2.4	CB1.1, CB2.1
2	Oligohaline	1–6	275	4.5	CB2.2, CB3.1
3	Mesohaline	5–11	362	6.6	CB3.2, CB3.3C
4	Mesohaline	9–15	909	10.2	CB4.1C, CB4.2C, CB4.3C, CB4.4
5	Mesohaline	12–17	1,475	10.5	CB5.1, CB5.1 W, CB5.2, CB5.3, CB5.4, CB5.4 W, CB5.5
6	Polyhaline	15–20	743	8.7	CB6.1, CB6.2, CB6.3, CB6.4
7	Polyhaline	18–23	1,521	8.9	CB7.1, CB7.1 N, CB7.1S, CB7.2, CB7.2E, CB7.3, CB7.3E, CB7.4 N, EE3.5
8	Polyhaline	20–25	412	7.7	CB7.4, CB8.1, CB8.1E
Bay		13–18	5,849	8.9	33 stations

Based on the mean annual cycle of the average salinity in each segment, 1998–2018.

From Chesapeake Bay Program (2004).

Computed as volume (Chesapeake Bay Program, 2004) divided by area.

Compared to its condition in the early 1900s, Chesapeake Bay is severely degraded due to eutrophication and the attendant declines in water clarity and dissolved oxygen (see Kemp et al., 2005, for a review). The reduction of water clarity, which may have been exacerbated by overharvesting of oysters and other filter feeders, has historically led to dramatic reductions in underwater grasses, an important benthic habitat. Over the past 30 years, Chesapeake Bay has experienced a modest improvement in overall water quality, as measured using metrics for chlorophyll a, dissolved oxygen, and water clarity/underwater grasses, consistent with expectations from reductions in riverine nitrogen loading (Zhang et al., 2018).

Studies of the carbonate system of the mainstem Bay are relatively recent. Observations indicate relatively poor acid buffering capacity (Cai et al., 2017); biological sources and sinks of DIC that vary with time, depth, and distance along the axis of the estuary (Brodeur et al., 2019; Friedman et al., 2020; Shadwick, Friedrichs, et al., 2019); and nonconservative TA behavior (Brodeur et al., 2019; Najjar et al., 2020). The only estimates of air‐water CO2 fluxes in the mainstem bay were based on an 8‐month data set in mesohaline and polyhaline waters (Friedman et al., 2020), a time series estimate at a single location in the polyhaline bay (Shadwick, Friedrichs, et al., 2019), and numerical models (Shen et al., 2019; St‐Laurent et al., 2020). These studies document considerable spatial and temporal variability in the flux, with the mainstem acting as both a source and sink of atmospheric CO2, depending on time and location.

Observations of pH, Temperature, Salinity, and Dissolved Oxygen

Observations of surface water pH, temperature, salinity, and dissolved oxygen are from the Chesapeake Bay Program (CBP), which conducts monthly (semimonthly in the warmer months) bay‐wide surveys. We selected 33 stations from the mainstem bay for the analysis (Figure 1 and Table 1), seeking good spatial coverage while avoiding overrepresentation of stations clustered together. Data from the beginning of the regular monitoring program in 1985 through 2018 were downloaded from the CBP's Water Quality Database (http://data.chesapeakebay.net). A close examination of the pH data led us to question some measurements in the early part of the record and limit our analysis to 1998–2018. The data were then subjected to quality control procedures to remove outliers. Supporting information Text S1 and Figure S1 provide details of the processing of pH, temperature, salinity, and oxygen data.

Total Alkalinity Mixing Model

TA is measured much less frequently than pH in Chesapeake Bay, and hence, we sought to model surface TA in the bay as a function of measured surface salinity and time. The relationship with salinity recognizes the nearly conservative nature of TA as well as the different TAs of the two main sources of water to the mainstem bay: Atlantic Ocean shelf waters (high‐salinity source) and the Susquehanna River (zero‐salinity source). The relationship with time recognizes the substantial seasonal and interannual variability of TA of the Susquehanna River. TA in the mainstem bay has increased over the past several decades at a rate that is approximately consistent with expectations from the increase in the riverine source and conservative mixing with a fixed marine endmember (see also Najjar et al., 2020). On the seasonal time scale, however, it is not clear to what extent river TA is felt in the mainstem bay. Hence, we used high‐quality measurements of TA from Shadwick, De Meo, and Friedman (2019), which are limited to 2016–2018, to develop and evaluate a TA mixing model made up of two components, one for S < 10, which accounts for seasonal and interannual variability in the riverine endmember, and one for S > 10, which accounts only for the long‐term increasing trend in TA in the mainstem bay. Supporting information Text S2 and Figure S2 provide the full details. Results, shown in Figure 2a, reveal a root‐mean‐square error of 129 mmol m−3 for S < 10 (n = 33) and 30 mmol m−3 for S > 10 (n = 135).

Figure 2

(a) TA predicted from the mixing model as a function of observed TA. (b) Predicted inorganic alkalinity as a function of observed TA. Observed TA is from Shadwick, De Meo, and Friedman (2019).

Organic Alkalinity Estimation

Organic alkalinity can lead to overestimation of pCO2 when pCO2 is computed from pH and TA. Hunt et al. (2011) found overestimation of 13–66% in North American rivers and Abril et al. (2015), in a synthesis of measurements from rivers in Africa, Amazonia, and Europe, found the median overestimation to be 23% but, importantly, showed that the overestimation increases dramatically with decreasing pH, increasing dissolved organic carbon concentration and decreasing TA. For example, the mean overestimation was 15% in rivers with pH > 7 (which have a mean pH of 7.9) and was 362% in rivers with pH < 7 (which have a mean pH of 6.3). Compared to rivers, the overestimation in coastal waters appears to be smaller, ranging from 10 to 47 μatm (roughly 2% to 10%) in the Baltic Sea (Kuliński et al., 2014), about 1% for typical conditions in Gulf of Mexico coastal waters (Yang et al., 2015), and −7% to 2% in the Scheldt Estuary (Frankignoulle & Borges, 2001).

Unfortunately, there are no readily available estimates of organic alkalinity in Chesapeake Bay. We estimated organic alkalinity by combining CBP pH measurements with simultaneous high‐quality measurements of DIC and TA (Shadwick, De Meo, & Friedman, 2019). The inorganic portion of TA was calculated using a MATLAB® version of the program CO2SYS (van Heuven et al., 2011) with inputs of pH, temperature, and salinity from the CBP and DIC from Shadwick, De Meo, and Friedman (2019). For this calculation, we used the carbonic acid equilibrium constants of Cai and Wang (1998), the equilibrium constant for bisulfate ion of Dickson (1990), and the ratio of total boron to salinity of Uppström (1974). Also, the phosphate and silicate concentrations were assumed to be zero as they have a minimal impact on the calculation. We chose the carbonic acid equilibrium constants of Cai and Wang (1998) because, as pointed out by Dinauer and Mucci (2017), these constants agree with the fresh water constants of Millero (1979) whereas other constants for low‐salinity waters (Millero, 2010; Millero et al., 2006) do not. The computed inorganic alkalinity is very close to the measured TA (Figure 2b) but is lower, on average, by 20 mmol m−3 and has a root‐mean‐square error of 30 mmol m−3. We therefore assumed the organic alkalinity in Chesapeake Bay is normally distributed with mean 20 mmol m−3 and ±standard deviation 30 mmol m−3. We considered making organic alkalinity a function of salinity but chose not to because the coefficient of determination is small (r 2 = 0.05).

Our estimate of organic alkalinity in Chesapeake Bay is within the range of literature values of organic alkalinity in rivers and coastal waters (see review in supporting information Text S5). More specifically, the value of 20 mmol m−3 is reasonably close to the value of 10 mmol m−3 found by Patsavas et al. (2015) for shelf waters of the eastern United States. An independent estimate of the freshwater endmember can be made using a model of organic alkalinity based on the concentrations of the hydrogen ion and dissolved organic carbon, which was developed from Baltic Sea data (Kuliński et al., 2014). Using a mean hydrogen ion concentration of 0.016 mmol m−3 in the tidal fresh region from our analysis and a Susquehanna River dissolved organic carbon concentration of 250 mmol m−3 (Fisher et al., 1998; Rochelle‐Newall & Fisher, 2002) leads to an organic alkalinity of 23 mmol m−3, which is very close to our estimate from DIC, TA, and pH data.

Calculation of Surface Water pCO2 and the Air‐Water CO2 Flux

Surface water pCO2 was computed from the observed surface pH, temperature, and salinity and the modeled TA minus the organic alkalinity (20 mmol m−3) at each station at the native temporal resolution of the CBP measurements (roughly monthly, section 2.1) using CO2SYS (van Heuven et al., 2011) with the assumptions described in section 2.4. Figure 3 shows that pCO2 increases with increasing temperature (T), decreasing salinity (S), decreasing pH, and increasing TA. Importantly, given the range of conditions in Chesapeake Bay, pCO2 is roughly an order of magnitude more sensitive to pH than to T, S, and TA. Furthermore, pCO2 is a highly nonlinear function of pH and a nearly perfect linear function of hydrogen ion concentration ([H+]), which becomes important when errors in the pH measurement are propagated in the pCO2 calculation.

Figure 3

Calculated surface water pCO2 as a function of (a) T, (b) S, (c) TA, (d) pH, and (e) [H+]. Note a factor of 10 larger vertical axis for pH and [H+]. Overlaid are histograms of the surface data. Both pH and [H+] are shown to emphasize differences in linearity. When T, S, TA, and pH (or [H+]) are allowed to vary, the other variables are held constant at their approximate mean values for the mainstem Chesapeake Bay surface waters (T = 18°C, S = 14, TA = 1.6 mol m−3, and pH = 8.1).

Monthly means of all quantities (pCO2, TA, temperature, salinity, and Δ[O2]) were computed for every station, gaps were filled, averages over each segment were computed, and finally, the air‐water CO2 flux was computed using a standard transfer‐velocity formulation (Wanninkhof, 2014) and winds from a reanalysis product. Unless specified otherwise, positive values indicate that the net CO2 flux is from water to air. All bay‐wide averages are computed by weighing segments by area. Full details are provided in supporting information Text S3.

Error Estimation

We assumed that the main uncertainties in monthly surface water pCO2 and the air‐water CO2 flux stem from errors in the pH measurement, the TA mixing model, the organic alkalinity estimation, the discrete sampling in time (once or twice per month), and the gas transfer velocity. A Monte Carlo analysis was conducted to propagate these errors to surface water pCO2 and the air‐water CO2 flux. Other errors, such as those associated with carbonate system equilibrium constants, are not included but are discussed in section 3.1.2.

pH and Alkalinity

For pH, the reported instrument error of the YSI and Hydrolab sondes is ΔpH = 0.2. This error was converted to an error in hydrogen ion concentration using Δ[H+] = 10–pH(10–ΔpH − 1), which is equal to 0.369[H+] for ΔpH = 0.2. The TA mixing model error was assumed to be the error in the fit of TA to S from Shadwick, De Meo, and Friedman (2019) high‐quality data set (Figure 2a, section 2.3). The organic alkalinity error was taken to be 20 ± 30 mmol m−3, as estimated in section 2.4 and Figure 2b. The errors in pH, the TA mixing model, and organic alkalinity were assumed fixed in time. The error distributions were assumed to be normal in [H+], TA, and organic alkalinity. Due to the strong nonlinearity of pCO2 as a function of pH (Figure 3a), sampling a normal distribution in pH results in a strongly skewed pCO2 distribution. The skewing is positive because the pCO2 versus pH curve is concave upward. Very little skewing occurs using [H+] and TA because pCO2 is essentially linear in these variables (Figure 2e). The resulting probability distributions were randomly sampled each time pCO2 was calculated, to create 1,000 realizations of pCO2 time series at each station in the bay. Occasionally, an unrealistic value of TA (negative) or [H+] (corresponding to pH less than 7 or greater than 9.69, which is outside of the observed pH range at the 33 stations used in the analysis, Figure 3d) was drawn from the distribution, in which case the value was discarded and another draw was made. The resulting time series, which matched the irregular temporal resolution of the Chesapeake Bay monitoring, were then averaged by month.

Discrete Sampling in Time

The monthly averages of pCO2 at each station are based on one or two pCO2 computations per month, so are not true monthly averages. To assess the associated sampling error, we used hourly surface pCO2 estimates at the mouth of the York River Estuary, which is in the polyhaline Chesapeake Bay (Shadwick, Friedrichs, et al., 2019). The pCO2 estimates were made by using direct measurements of pH, T, and S and a model of TA as a function of S. The data span the time periods December 2016 to March 2017 and January 2017 to May 2018, a total of nine complete months. The time periods (winter and spring) and single location of this data set are not ideal but provide the opportunity to make a first estimate of pCO2 sampling error in Chesapeake Bay and, as far as we are aware, in any estuary. While Shadwick, Friedrichs, et al. (2019) first computed diurnal averages of pH, T, and S before computing pCO2, we retain the original hourly resolution in order to capture diurnal and tidal influences. Our computation also differs from that of Shadwick, Friedrichs, et al. (2019) in that we use the equilibrium constants of Cai and Wang (1998) instead of those of Mehrbach et al. (1973) refit by Dickson and Millero (1987).

To compute the sampling error, we first created a probability distribution at hourly resolution of the times of day that pH was measured in our version of the CBP data set. This probability distribution was then used to randomly sample the hourly pCO2 data set, assuming no preference on the day of the month. The difference between sampled pCO2 and the true mean pCO2 for the month that was sampled was determined to be the sampling error (vertical deviation from the 1:1 line in Figure 4). A sampling error distribution was created by repeating the sampling process 1,000 times. This error distribution has a mean of −9.6 μatm and a standard deviation of 49 μatm. For each of the 1,000 pCO2 realizations at the 33 stations in the mainstem Bay (section 2.5), the pCO2 sampling error was added by drawing from the sampling error distribution. If a draw resulted in a negative pCO2, the pCO2 was set to zero.

Figure 4

Ensemble of instantaneous surface water pCO2 draws as a function of true monthly mean surface water pCO2 based on the high‐temporal resolution data set of Shadwick, Friedrichs, et al. (2019).

Gas Transfer Velocity

Wanninkhof (2014) attempted to quantify the uncertainty in the gas transfer velocity due to the combined influences of uncertainties in the leading coefficient a, the Schmidt number, winds, high‐wind speed effects, and chemical enhancement. The result was an overall uncertainty of 20%, which is applicable to global‐scale estimates of air‐sea CO2 exchange. This is not a rigorous estimate, and it is unclear if 20% is one standard error or half the 95% confidence range. Quantifying the error at the local scale in a rigorous manner is not possible at this time. To be conservative, we used a standard error estimate of 30%. Specifically, every time the gas transfer velocity was computed, which is once per month per segment, we used a Monte Carlo approach and randomly added an error chosen from a normal distribution with a mean of 0 and standard deviation of 30%.

Error Reporting Conventions

We use the following convention to report the pCO2 and CO2 flux calculation results. The results of a single pCO2 and CO2 flux calculation (i.e., without random errors in pH, the alkalinity mixing model, organic alkalinity, temporal sampling, and the gas transfer velocity) are referred to as point estimates. The results of Monte Carlo calculations are reported as best estimates (means) with 95% confidence intervals (2.5th and 97.5th percentiles) or standard errors (standard deviations) of each corresponding ensemble. Unless otherwise stated, all results reported without uncertainties are point estimates, and those reported with a confidence range or standard error are best estimates.

Results

pCO2 Variability

Temporal Variability

Surface water pCO2 is highly variable in time, as time series (Figure 5) and basic statistics (Table 2) for each segment reveal. Segments 1 and 2 (tidal fresh and oligohaline regions) stand out as particularly variable in surface water pCO2, with the highest coefficients of variation, 66% and 65%, respectively, compared to 34–57% in the rest of the bay. Similarly, Segments 1 and 2 have the most prominent annual cycles, making up 56% and 69%, respectively, of the variability of the full time series, compared to 14–33% in the mesohaline and northern polyhaline bay (segments 3–6) and <5% in the southern polyhaline bay (Segments 7–8).

Figure 5

Temporal variability of the best estimates of the surface water pCO2 computed for eight segments (Figure 1) of the mainstem Chesapeake Bay.

Table 2

Statistics of the Best Estimates of Monthly Average Surface Water pCO2 in Each Segment and the Bay‐Wide Average

Segment	Mean ± 1 standard deviation (μatm)	Coefficient of variation	Mean annual cycle
			Variance fraction	Month of minimum	Month of maximum
1	803 ± 534	0.66	0.56	February	September
2	1,139 ± 740	0.65	0.69	February	July
3	576 ± 309	0.54	0.22	March	September
4	378 ± 192	0.51	0.28	March	October
5	345 ± 145	0.42	0.33	May	September
6	373 ± 214	0.57	0.14	May	October
7	449 ± 196	0.44	0.04	May	September
8	527 ± 244	0.46	0.01	January	August
Bay	457 ± 154	0.34	0.25	March	September

Figure 6 shows that mean annual cycle is well defined in Segments 1 and 2 (tidal fresh and oligohaline region) and Segments 4–7 (central mesohaline to central polyhaline region), with the saltier segments lagging behind the fresher segments by 0–3 months. From the fresher segments to the saltier segments, minima shift from winter to typically spring and maxima shift less, from summer to sometimes early fall (Table 2).

Figure 6

Mean annual cycles of the best estimates of surface water pCO2 and oxygen supersaturation (Δ[O2]) in eight segments (Figure 1) of the mainstem Chesapeake Bay. Note the change in vertical axes among the panels. Dashed lines show atmospheric pCO2 (blue) and oxygen saturation (red). Gray shading shows 95% confidence range on pCO2.

Biological control on the temporal variability in surface water pCO2 is suggested by the anticorrelation between pCO2 and Δ[O2], which is clearly seen in mean annual cycles of these two quantities (Figure 6). The extrema in the two annual cycles line up well. The lag of the lower bay with respect to the upper bay that is seen in pCO2 is also seen in Δ[O2]. Table 3 shows that Δ[O2] accounts for 47% to 90% of the variability in the annual cycle of pCO2 in six of the eight segments. Table 3 also shows >95% significance in the relationship between the two variables for the full time series in all segments, with 25% to 43% of the variability accounted for in Segments 1–5. Even with the mean annual cycle removed, all segments show correlations of greater than 95% confidence between the two variables. Only for annual averages are correlations generally poor.

Table 3

Coefficient of Determination (r 2) of Surface Water pCO2 (Best Estimates) Versus Oxygen Supersaturation and Salinity for Different Time Frames

Segment	Oxygen supersaturation				Salinity
	Full time series	Annual averages	Anomalies	Mean annual cycle	Full time series	Annual averages	Anomalies	Mean annual cycle
1	0.25	0.33	0.18	0.47	0.11	0.36	0.08
2	0.35	0.21	0.14	0.81	0.05	0.45	0.06
3	0.37	0.23	0.42
4	0.43		0.35	0.74	0.15	0.21	0.06	0.64
5	0.31		0.16	0.80	0.39	0.69	0.28	0.75
6	0.09		0.03	0.90	0.20	0.22	0.11	0.87
7	0.03		0.02	0.55	0.20	0.32	0.20	0.47
8	0.04		0.03		0.10		0.12
Bay	0.18		0.09	0.66	0.25	0.48	0.23

Note. Anomalies are the full time series with the mean annual cycle removed. Values are shown when there is >95% confidence that the slope of the relationship is different from zero; under these conditions, all correlation coefficients (r) are negative with oxygen supersaturation and positive with salinity.

Streamflow has the potential to profoundly impact estuarine pCO2 through its influence on salinity, stratification, and the delivery of inorganic nutrients, organic matter, exogenous DIC, and sediment. To investigate the potential influence of streamflow on pCO2, correlation statistics are presented with salinity (Table 3). We prefer using salinity instead of streamflow in the correlation analysis because it takes time, sometimes several months, for a change in flow and other characteristics of the Susquehanna River to be felt throughout the mainstem bay (Shen & Wang, 2007). There is a close correspondence between pCO2 and salinity for the four time scales considered in Table 3. All of the significant correlations are positive. Correlations are particularly strong for annual averages, where 21–69% of the variability in pCO2 is accounted for by salinity in six of the eight segments. Correlations are also strong for the annual cycle in Segments 4–7, where salinity accounts for 47–87% of the variability in pCO2. A positive correlation between salinity and pCO2 for constant temperature, DIC, and TA is expected through the salinity dependence of the equilibrium constants, but the sensitivity found here is higher. For example, for the bay‐wide average correlations (last row in Table 3), the slopes are equivalent to 6–7% per salinity unit, while the expectation from the carbonate equilibria is 4% per salinity unit with T, S, and DIC and TA at their bay‐wide average values of 15.4°C, 15.3, 1,462 mmol m−3, and 1,537 mmol m−3, respectively. The positive pCO2‐salinity correlation indicates that periods of high flow are associated with lower pCO2. This is possibly a result of riverine nutrient delivery, which would enhance photosynthesis and decrease pCO2. In a net autotrophic system like Chesapeake Bay (Kemp et al., 1997), this effect may have more of an influence than the simultaneous delivery of sediment and organic matter, which would tend to push the estuary toward heterotrophy via light limitation by sediment and enhanced respiration due to organic matter.

Spatial Variability

Spatial variability in pCO2 is considerable, with the best estimates for the long‐term averages varying among segments by a factor of 3 (Table 2 and Figure 7): from a low of 345 (336, 352) μatm in Segment 5 (mesohaline) to a high of 1,139 (1,098, 1,177) μatm in Segment 2 (oligohaline). Similar to the temporal variability of pCO2, the spatial variability of pCO2 appears to be biologically driven, as there is a strong anticorrelation with Δ[O2] across the eight segments (r 2 = 0.91, p < 0.001, Figure 7).

Figure 7

Long‐term averages of the best estimates of surface water pCO2 and oxygen supersaturation (Δ[O2]) in eight segments (Figure 1) of the mainstem Chesapeake Bay as a function of latitude. Dashed lines show atmospheric pCO2 (blue) and oxygen saturation (red). Error bars show 95% confidence range on pCO2.

Degree of Supersaturation

Figures 5, 6, 7 show that bay surface waters vary in their degree of saturation with respect to atmospheric CO2. The bay‐wide, long‐term average surface water pCO2 is 457 (452, 462) μatm, which is 17% greater than the corresponding atmospheric value, 381 μatm. Segments vary from being, on average, highly supersaturated (199% in Segment 2) to slightly undersaturated (−10%, Segment 5). Undersaturation occurs in all segments, but the fraction of time this occurs varies greatly, as seen in Figure 5, from only 9% in Segment 2 to 66% in Segment 6. The month of most frequent undersaturation by segment is always in winter or spring (January–June). As can be seen in Figures 5 and 6, the seasonality of the degree of saturation is almost completely controlled by surface water pCO2, which has a range (maximum minus minimum) in its mean annual cycle that varies from 97 μatm in Segment 8 to 1,628 μatm in Segment 2, much greater than the range in the mean annual cycle of atmospheric pCO2 (15–17 μatm among the eight segments).

Air‐Water Flux

The air‐water CO2 flux is the product of the pCO2 gradient across the air‐water interface, the CO2 gas transfer velocity, and the CO2 gas solubility (Equation S3). We have already noted that variability in the pCO2 gradient is driven almost completely by surface water pCO2 (section 3.1.3). Temperature has a large impact on CO2 gas solubility () and, through the Schmidt number, on the gas transfer velocity (k ). However, temperature effects on the product k nearly cancel, and so it is convenient to define this product as the gas transfer coefficient, k , whose variability is mainly due to wind speed (Etcheto & Merlivat, 1988), with some minor variability due to salinity and total air pressure.

In the mainstem bay, k has considerable temporal (Figures S3a and S3b) and spatial (Figure S3c) variability, both of which, as expected from the formulation of the gas transfer velocity adopted here (Equation S4), are dominated by the variability in the square of the wind speed. Temporal variability in the bay‐averaged k is largely seasonal, with the mean annual cycle accounting for 66% of the variability of the full time series and characterized by a maximum in January that is 2.5 times higher than the minimum, which occurs in July (Figure S3b). There is also a modest but significant (p = 0.02) declining trend in the bay‐averaged k anomaly of −9.04 × 10−6 mol m−2 month−1 μatm−1 yr−1, which corresponds to a decline of 11% from 1998 to 2018 (given the mean k of 0.00174 mol m−2 month−1 μatm−1). Temporal variability in individual segments is similar to the bay‐wide average, with strong annual cycles that account for 58% to 76% of the variability of the full time series. Spatial variability of long‐term averages shows a slight decline from Segments 1 to 4, an increase of 47% from Segments 4 to 5, and a modest decline from Segments 5 to 8.

Averaged over the bay, the air‐water CO2 flux is highly variable, with a prominent annual cycle, and fewer periods of CO2 outgassing toward the end of the study period (Figure 8a). To a large extent, the variability of the bay‐averaged air‐water CO2 flux mimics the variability of the bay‐averaged air‐water pCO2 gradient, with both having strong mean annual cycles with spring minima and fall maxima (Table 2 and Figure 8b). Spatial variability in the CO2 flux, with high CO2 outgassing in the upper bay and weak‐to‐moderate CO2 outgassing in the lower bay (Figure 8c), is also largely driven by spatial variability in the air‐water pCO2 gradient (Figure 7).

Figure 8

Best estimates of the net flux (with 95% confidence intervals) of CO2 from water to air for Chesapeake Bay. Bay‐wide average (a) monthly (blue) and annual (red) time series and (b) mean annual cycle. (c) Long‐term averages in flux per unit area (blue) and integrated flux (red). The blue x's are the mean gas transfer coefficient multiplied by the mean pCO2 gradient across the air‐water interface.

On average the bay outgasses CO2 at rate of 1.2 (1.1, 1.4) mol m−2 yr−1. The mean flux is only slightly less than the product of the mean gas transfer coefficient (0.0209 mol m−2 yr−1 μatm−1) and the mean air‐water pCO2 gradient (66 μatm), which is equal to 1.4 mol m−2 yr−1. Thus, covariations between the gas transfer coefficient and the air‐water pCO2 gradient are modest, due to the fact that most of the variability in k is seasonal and 90° out of phase with the bay‐wide average air‐water pCO2 gradient (e.g., the maximum in k is in winter and the maximum of the gradient is in fall). In fact, the product of the mean k and mean air‐water pCO2 gradient departs substantially from the mean flux only in Segments 1 and 2 (Figure 8c). In these segments, this product exceeds the mean flux because the mean flux accounts for the negative covariation between k and the air‐water pCO2 gradient, which results from the mean annual cycles of these two quantities being nearly 180° out of phase (e.g., compare the air‐water pCO2 gradient in Figure 6a with k in Figure S3b).

The air‐water CO2 flux integrated over each segment has a substantially different pattern with latitude than the flux per unit area (Figure 8c) because of the large range in the areas of the eight segments (Table 1); for example, the largest segment (7) is 10 times the area of the smallest segment (1). The bay‐wide integrated flux is 84.9 Gg yr−1, with Segments 1 + 2 + 3 outgassing 58.8 Gg yr−1, Segments 4 + 5 ingassing 16.0 Gg yr−1, and Segments 7 + 8 outgassing 41.1 Gg yr−1.

The bay‐averaged annual air‐water CO2 fluxes vary considerably over the study period, from −0.8 (−1.3, −0.4) mol m−2 yr−1 in 2006 to 7.2 (6.1, 8.1) mol m−2 yr−1 in 1999 (Figure 8a). Closer examination of the annual fluxes reveals that none of the annual fluxes during the first 5 years (1998–2002) were exceeded in the following years. These unusually high annual fluxes, with an average (±1 standard deviation) of 4.1 ± 1.9 mol m−2 yr−1, are a stark contrast to the rest of the record, which has an average flux of 0.3 ± 1.0 mol m−2 yr−1. These high fluxes are related to relative high salinity (16.6 ± 1.9) compared to the rest of the record (15.3 ± 1.7), which is due to the relatively low streamflow (Figure S2a). As discussed in section 3.1.1, a net autotrophic system such as Chesapeake Bay might be expected to be pushed toward more heterotrophy (and hence higher CO2 outgassing) under low‐flow years.

Uncertainty Analysis

We accounted for errors in the pH measurement, the alkalinity mixing model, the organic alkalinity estimation, temporal sampling, and the gas transfer velocity and now compare them for the calculation of the long‐term, bay‐wide averages of the surface water pCO2 and air‐water CO2 flux (Figure 9). Six cases are considered: all errors and each of the five errors alone. The best estimate is always slightly higher than the point estimate as a result of weak nonlinearities in the surface water pCO2 calculation. The biggest error in the surface water pCO2 calculation is the pH measurement error (Figure 9a). The same is true of the air‐water CO2 flux calculation, with errors, the gas transfer velocity, temporal sampling, the alkalinity mixing model, and the organic alkalinity estimation having 95% confidence ranges in surface water pCO2 that are 72%, 27%, 15%, and 5% of the same range due to the pH measurement error (Figure 9b). Errors may be large for an individual month for an individual segment, but they decrease dramatically when averaged over time and space. Figure 10a shows an example that starts with a fairly large error (95% confidence range) in the surface pCO2 for November 2009 in Segment 7 (290 μatm) and gets progressively smaller for the 2009 average in that segment (117 μatm), the long‐term average in that Segment (20 μatm), and then finally for the long‐term, bay‐wide average (10 μatm). Similar reductions are seen for the air‐water CO2 flux (Figure 10b). The dramatic error reduction due to averaging emphasizes one of the advantages of water quality monitoring in Chesapeake Bay, which is characterized by a long and regular record with excellent spatial and temporal coverage.

Figure 9

Best estimate and 95% confidence ranges with the point estimate removed for the long‐term, bay‐wide averages of the (a) surface water pCO2 and (b) air‐water CO2 flux. Six cases are considered: all errors, the pH measurement error, the alkalinity mixing model error, the organic alkalinity estimation error, the temporal sampling error, and the gas transfer velocity error.

Figure 10

Example of the effect of temporal and spatial averaging on the error (95% confidence range) in (a) surface water pCO2 and (b) the air‐water CO2 flux. The bars in each panel represent, in order, segment 7 November 2009 monthly average, segment 7 2009 annual average, segment 7 long‐term average, and the bay‐wide long‐term average.

Discussion

Other Uncertainties

Uncertainties in Carbonate System Equilibria

In addition to uncertainties in pH, the TA mixing model, organic alkalinity, temporal sampling, and the gas transfer velocity, there are uncertainties in the equations of the carbonate system equilibria, which are used to compute pCO2 from pH and TA. Such uncertainties stem from the assumption of zero phosphate and silicate concentration and uncertainties in carbonate system equilibrium constants. Each of these is considered now.

Phosphate and silicate concentrations in surface waters of Chesapeake Bay can be as high as 1.5 and 100 μmol kg−1, respectively (Malone et al., 1996). For the mean conditions described above (Figure 3, section 2.5), changing phosphate and silicate concentrations to these concentration decreases surface water pCO2 by only 0.4 and 0.6 μatm, respectively. Hence, the assumption of zero concentration for these nutrients is reasonable.

Random errors in the carbonic acid equilibrium constants also translate to small errors in pCO2, about 2 μatm, when computed from pH and TA (Emerson & Hedges, 2008, page 115). However, systematic errors may be substantial, particularly in low‐salinity waters (Dinauer & Mucci, 2017). In addition to the carbonic acid constants of Cai and Wang (1998), which were used here and are applicable for S = 0–40, the CO2SYS program (van Heuven et al., 2011) has two other sets of equilibrium constants applicable for estuarine waters: Millero (2010) for S = 0–40 and Millero et al. (2006) for S = 1–50. We compared best estimates of the air‐water CO2 flux and surface water pCO2 computed with these three sets of constants as well as the constants of Mehrbach et al. (1973) as refit by Dickson and Millero (1987), which are applicable for 20 < S < 40, are widely believed to be the most accurate for oceanic waters (e.g., Emerson & Hedges, 2008), and gave good results when comparing pCO2 computed from pH and TA with measured pCO2 in estuarine conditions (Frankignoulle & Borges, 2001). The bay‐wide, long‐term average CO2 outgassing using the constants of Millero (2010), Millero et al. (2006), and Mehrbach et al. (1973) was 47%, 46%, and 23% higher, respectively, than that using the constants of Cai and Wang (1998). Surface water pCO2 computed using the constants of Millero (2010) was higher than that computed using the constants of Cai and Wang (1998) by an average of 305 μatm in Segment 1 and decreasing dramatically and monotonically to 5 μatm in Segment 8, consistent with the findings of Dinauer and Mucci (2017). Differences between the constants of Millero et al. (2006) and Millero (2010) were everywhere small, with average differences in surface water pCO2 of 0–2 μatm in all segments. Surface water pCO2 computed with the constants of Mehrbach et al. (1973) was lower in Segments 1 and 2 (114 and 20 μatm, respectively) and slightly higher in Segments 3–8 (9–23 μatm) than that computed with the constants of Cai and Wang (1998).

Additional Uncertainties in the Gas Transfer Velocity

The gas transfer velocity parameterization used here (Wanninkhof, 2014) was designed primarily for use in the open ocean. The main difference in the gas transfer velocity between estuaries and the open ocean is that a significant fraction of the turbulence that regulates gas exchange in estuaries may originate from shear at the sea floor resulting from tidal currents (Ho et al., 2016). However, we show in supporting information Text S4 that it is unlikely that tidal currents have a significant impact on the gas transfer velocity in Chesapeake Bay, which is relatively deep and has relatively weak tidal currents. There are other gas transfer processes distinct to estuaries, including fetch and turbidity. However, the impact of these factors is highly uncertain because parameterizations of these processes are largely untested (Borges & Abril, 2011).

Comparison With Other Estuarine Air‐Water CO2 Flux Studies

We compared our air‐water CO2 flux results with three recent studies in the mainstem bay. Shen, Testa, Li, et al. (2019) employed a biogeochemical model to estimate the air‐water CO2 flux in the mainstem bay for the year 2016 and found a net uptake of 2.2 mol C m−2 yr−1 or 130 Gg C yr−1 over an area of 5,050 km2, somewhat smaller than the area used here. Shen, Testa, Ni, et al. (2019) extended the model simulation to the 1985–2016 period and found that the average ±1 standard deviation of the annual mainstem bay fluxes was −39 ± 112 Gg C yr−1, indicating uptake, on average, but with large interannual variability. Friedman et al. (2020) estimated the air‐water CO2 flux from surface water observations of pCO2 during four cruises from fall 2016 to summer 2017 in the mesohaline and polyhaline bay to find a net uptake of 0.38 mol m−2 yr−1, which translates to 23 Gg C yr−1 using the sum of areas 4–8 in Table 1. Shen, Testa, Ni, et al. (2019) found relatively high rates of CO2 outgassing in the upper bay (Segments 1–3) during 2016: 8.8 mol C m−2 yr−1, which agrees well with our estimate of 10.0 (7.5, 12.8) mol C m−2 yr−1 for the same region and year. In the middle bay (Segments 4 and 5), Shen, Testa, Li, et al. (2019) and Friedman et al. (2020) found an uptake of approximately 3 mol C m−2 yr−1 during 2016 and during a time period from fall 2016 to summer 2017, respectively. These estimates contrast with our estimate of a weak CO2 outgassing at a rate of 0.7 (0.0, 1.5) mol C m−2 yr−1 in 2016 and near‐neutral flux of −0.2 (−0.9, 0.6) mol C m−2 yr−1 from fall 2016 through summer 2017 (negative flux indicates uptake). In the lower bay (segments 6–8) in 2016 we estimated near neutral condition of–0.1 (−0.8, 0.6) mol C m−2 yr−1, which is in agreement with the findings of Shen, Testa, Li, et al. (2019) who also found neutral conditions in 2016. In the same geographic region but during a period from fall 2016 to summer 2017, Friedman et al. (2020) found and uptake of about 0.5 mol C m−2 yr−1, which is slightly outside our flux range of between −0.3 and 1.1 mol C m−2 yr−1.

Although Chesapeake Bay is typical of other estuaries in that it is a source of CO2 to the atmosphere, it is a very weak source. In the global synthesis of Chen et al. (2013), the mean area‐weighted estuarine CO2 outgassing rate is 7.7 mol C m−2 yr−1, a factor of 7 greater than the mean outgassing of Chesapeake Bay found in our study (1.1–1.4 mol C m−2 yr−1). The bay is also is a weak source compared to estuaries of eastern North America, which outgas CO2 at a rate of 9 ± 4 (mean ± 2 standard errors) mol C m−2 yr−1. Chen et al. (2013) found the mean CO2 outgassing for S < 2 to be 39.0 mol C m−2 yr−1 and for 2 > S > 25 to be 17.5 mol C m−2 yr−1, which can be compared to Chesapeake Bay's mean outgassing rates for Segments 1–2 and 3–8: 8.6 and 0.5 mol C m−2 yr−1, respectively. Hence, Chesapeake Bay follows the global outgassing pattern with salinity, albeit at much lower rates. The relative weakness of Chesapeake Bay's outgassing is consistent with it being net autotrophic, while most of the world's estuaries are heterotrophic (Borges & Abril, 2011). Note that, due to lateral advection, a net autotrophic system need not be a sink of atmospheric CO2.

Inorganic Carbon Budget of the Mainstem Chesapeake Bay

The results of this study in combination with studies on the riverine input of DIC and NEP allow an inorganic carbon budget to be developed for the mainstem bay. If averages over a year or more are considered, then it is reasonable to assume that the mainstem bay is in steady state, with riverine input (R) and lateral transport (L) of DIC balanced by losses due to NEP and air‐water exchange of CO2 (f): where the dimensions of the equation are mass per unit time. Here we have assumed that sources and sinks of calcium carbonate are relatively small, which was suggested by the near‐neutral alkalinity budget of the mainstem bay developed by Brodeur et al. (2019). Lateral transport represents advection and turbulent diffusion of DIC at the boundaries of the mainstem bay with Atlantic Ocean waters at the mouth of the bay, major tributaries of the mainstem bay (e.g., the Potomac, James, York, and Rappahannock River estuaries), and fringing marshes. Lateral transport is difficult to measure and is computed here as a residual term of Equation 1 after estimates of R, NEP, and f are made.

Riverine input of DIC to the Chesapeake Bay mainstem from the Susquehanna River was estimated by Stets and Striegl (2012) to be 0.27 Tg C yr−1. We update this flux here using the approach for generating the river alkalinity time series (supporting information Text S2). A time series of daily DIC fluxes in the Susquehanna River at Conowingo, Maryland (USGS gauge number 01578310; USGS, 2019) was constructed for 1 October 1967 to 30 September 2016 by first computing DIC concentrations from pH and TA, as in Stets and Striegl (2012), and then applying a statistical model (Hirsch et al., 2010); effective monthly average DIC concentrations were computed as the monthly fluxes divided by the monthly streamflow. To gap fill effective monthly average concentrations for October 2016 to December 2018, DIC was modeled as a function of streamflow (same as for river alkalinity in supporting information Text S2). DIC fluxes at monthly resolution for the period of interest here, 1998–2018, were calculated as a product of effective monthly concentrations and monthly flows and then averaged annually. The mean annual flux (±1 standard deviation, taken here to represent the uncertainty on the long‐term annual mean) is equal to 0.42 ± 0.13 Tg C yr−1. This is higher than the flux estimated by Stets and Striegl (2012) because of the different time period covered by that study.

There are several estimates of NEP in the mainstem bay. Kemp et al. (1997) employed five independent methods based on incubations and the budgets of carbon, nitrogen, and phosphorus, finding the mean of these approaches (± its standard error) to be 4.2 ± 0.6 mol C m−2 yr−1 or 0.29 ± 0.04 Tg C yr−1 using the area in Table 1. Much of the data that went into these estimates are from the late 1970s to early 1990s. Feng et al. (2015) applied a nitrogen‐based biogeochemical model over the period 2001–2005 to find annual NEP (±1 standard deviation of the five annual averages) to be 74 ± 23 × 109 g N yr−1, which translates to 0.46 ± 0.14 Tg C yr−1 using a C:N ratio of 106:17 (Hedges et al., 2002). Finally, Friedman et al. (2020) estimated NEP to be −0.5 mol C m−2 yr−1 by constructing a DIC budget in the mesohaline and polyhaline mainstem bay from fall 2016 to the summer of 2017. The range in these estimates is not surprising given the different approaches, time periods, and spatial extent of the studies. Given the limited spatial and temporal extent of the Friedman et al. (2020) estimate, we do not consider it further here. Rather, we take the average of the estimates of Kemp et al. (1997) and Feng et al. (2015) to arrive at a mean NEP (±1 standard error) of 0.38 ± 0.15 Tg C yr−1, where we assumed that the standard deviation of the Feng et al. (2015) estimate represents the standard error and that the errors from the two studies are independent.

Combining our estimated net flux of CO2 from the mainstem bay to the atmosphere of 0.085 ± 0.005 Tg C yr−1 with the estimates of riverine input and NEP described above and Equation 1 leads to a net lateral input of DIC to the mainstem bay of 0.04 ± 0.20 Tg C yr−1. We have computed the uncertainty in L by propagating the above uncertainties in R, NEP, and f, assuming the uncertainties to be independent of each other. The positive value of L indicates that there is likely a source, albeit weak, of DIC to the mainstem bay from lateral advection. This source, for example, could be due to the residual circulation at the mouth of the bay, which transports DIC‐rich water into the bay at depth and DIC‐poor water out of the bay near the surface. Our findings can be compared with those of Brodeur et al. (2019), who applied a 1‐D, steady‐state mixing model to 12 monthly surveys of the mainstem bay and found that the bay exports DIC at rate of 0.29 ± 0.06 Tg C yr−1 at its mouth. Our findings are not necessarily inconsistent with those of Brodeur et al. (2019) because of the different time periods considered as well as the other locations where DIC can enter the mainstem bay (e.g., from the tributaries and fringing marshes).

A local view of the carbon balance is afforded by comparisons of NEP with the air‐water CO2 flux in three main portions of the bay that correspond fairly closely to those used by Kemp et al. (1997) in their NEP analysis based on incubations: the upper, middle, and lower bay, which encompass Segments 1–3, 4–5, and 6–8, respectively (Figure 1). In the upper bay, Kemp et al. (1997) found NEP to be −7.3 mol C m−2 yr−1, indicating net heterotrophic conditions or net production of CO2. For our study period, we find that the mean net CO2 flux from the estuary to the atmosphere in the upper bay is 6.2 (5.8, 6.6) mol C m−2 yr−1, which nearly balances Kemp et al. 's (1997) NEP in the upper bay, indicating that net lateral fluxes of DIC into or out of this region are small. In the middle bay, there is rough agreement between Kemp et al. 's (1997) NEP of 0.0 mol C m−2 yr−1 and our weak CO2 flux of −0.5 (−0.7, −0.3) mol C m−2 yr−1, both of which indicate near‐neutral conditions. The lower bay is where the two fluxes diverge the most, with Kemp et al. 's (1997) NEP equal to 7.7 mol C m−2 yr−1 and our CO2 outgassing of 1.3 (1.1, 1.5) mol C m−2 yr−1, which implies, at steady state, a net lateral transport into the lower bay of about 9 mol C m−2 yr−1.

In summary, we find that the air‐water CO2 flux is a modest component of the inorganic carbon budget for the mainstem bay as a whole but can be locally significant. CO2 outgassing accounts for only about a fifth of the net DIC sinks in the mainstem bay (CO2 outgassing + NEP). Locally, the CO2 outgassing flux may either be the main way DIC produced from NEP is removed from the system (e.g., in the upper bay) or a small fraction of the net lateral influx (e.g., in the lower bay).

Summary and Conclusions

An analysis of surface water pCO2 and the air‐water CO2 flux has been conducted for eight segments of the mainstem Chesapeake Bay at monthly resolution for a 21‐year period (1998–2018) by exploiting pH, TA, temperature, and salinity measurements from a water quality monitoring program. To our knowledge, this is the first attempt to simultaneously propagate errors in multiple factors, including pH, TA, organic alkalinity, temporal sampling, and the gas transfer velocity. The main findings are as follows:

The largest error results from the measurement of surface water pH, closely followed by the error in the gas transfer velocity. Errors in the alkalinity mixing model and temporal sampling were secondary, but substantial, while the error in organic alkalinity was estimated to be very small.

On average the bay is a weak source of CO2 to the atmosphere, with a net flux of 1.2 (1.1, 1.4) mol m−2 yr−1 or 85 (75, 95) Gg yr−1.

Surface water pCO2 is highly variable in time, particularly in the tidal fresh and oligohaline bay, where there is a coefficient of variation as high as 66% and a fraction of the variability accounted for by the mean annual cycle as high as 69%.

There is a shift in the annual cycle of surface water pCO2 from the fresher segments to the saltier segments of the bay, with minima shifting from winter to typically spring and maxima shifting less, from summer to sometimes early fall.

Biological control on the seasonal, interannual, and spatial variability in surface water pCO2 is indicated by the strong anticorrelation between pCO2 and oxygen supersaturation.

pCO2 and salinity are positively and significantly correlated over seasonal and interannual timescales, even after correcting for the salinity dependence of the carbonate system equilibrium constants, suggesting that periods of high streamflow lead to more autotrophy.

Although historic pH and alkalinity measurements from water quality monitoring programs are generally of lower quality than those now made routinely by the marine chemistry research community, the shear abundance of them allows random measurement errors to be reduced to the point where meaningful insights can be gained about spatial and temporal variability in important derived quantities, such as surface water pCO2 and the air‐water CO2 flux. We demonstrated that high‐quality and high‐resolution carbonate system measurements of limited spatial and temporal coverage can be paired successfully with the abundant historical measurements of pH and salinity to facilitate bay‐wide analysis of carbonate system over multiple decades. Our findings suggest that the abundance of historical pH measurements in estuaries around the globe should be mined in order to constrain the large spatial and temporal variability of the CO2 exchange between estuaries and the atmosphere, leading to an improved understanding of the role of estuaries in the global carbon cycle.

Supporting Information S1

Click here for additional data file.