Redox reactions and weak buffering capacity lead to acidification in the Chesapeake Bay

The combined effects of anthropogenic and biological CO2 inputs may lead to more rapid acidification in coastal waters compared to the open ocean. It is less clear, however, how redox reactions would contribute to acidification. Here we report estuarine acidification dynamics based on oxygen, hydrogen sulfide (H2S), pH, dissolved inorganic carbon and total alkalinity data from the Chesapeake Bay, where anthropogenic nutrient inputs have led to eutrophication, hypoxia and anoxia, and low pH. We show that a pH minimum occurs in mid-depths where acids are generated as a result of H2S oxidation in waters mixed upward from the anoxic depths. Our analyses also suggest a large synergistic effect from river–ocean mixing, global and local atmospheric CO2 uptake, and CO2 and acid production from respiration and other redox reactions. Together they lead to a poor acid buffering capacity, severe acidification and increased carbonate mineral dissolution in the USA’s largest estuary.

A nthropogenic carbon dioxide (CO 2 ) has increased more rapidly in the atmosphere since the Industrial Revolution than natural CO 2 increase in any period of the last 800,000 years 1, 2 ; consequently, it has been known that the uptake of CO 2 by the ocean has altered surface seawater acidbased chemistry lowering pH by about 0.1 unit and calcium carbonate saturation state by roughly 0.5. This process, known popularly as ocean acidification (OA) for over a decade, will continue to decrease seawater pH by about 0.3 units by the end of the century 3,4 . It is likely that OA will cause detrimental effects on the health of marine organisms and ecosystems and alter the associated biogeochemical processes [5][6][7] .
Recent research indicates that eutrophication can exacerbate OA, where respiratory processes contribute a far greater acidification in the coastal oceans relative to the open ocean [8][9][10][11][12][13] . Coastal eutrophication occurs with increased inputs of nutrients from the application of chemical fertilizers, discharges of human and animal wastes, and atmospheric NO x inputs from fossil fuel burning, which have fueled large algal blooms in many coastal water bodies, especially those near population centers 14 . It is well known that decomposition of algal organic matter from highly productive surface water leads to the development of seasonally low oxygen (hypoxic) or even zero oxygen (anoxic) bottom waters in many coastal water bodies in the world 15,16 . However the coupling between redox and acid-base chemistry has not been explored extensively in seasonally anoxic and partially mixed estuaries nor in permanently anoxic deep basins although redox chemistry and pH have been reported before in the latter 13,[17][18][19][20][21][22] . Specifically, it is not known how subsurface water pH dynamics are influenced by anaerobic respiration and the oxidation of reduced chemical species (notably H 2 S) in seasonally low oxygen (O 2 ) estuaries around the world let alone the interaction of these processes with the anthropogenic CO 2 induced OA.
The Chesapeake Bay is the largest estuary in the United States with a well-documented history of eutrophication over the past half century [23][24][25] . A recent report demonstrates that some regions of the bay have suffered a long-term pH decline related to eutrophication 26 . However, few process studies have examined the CO 2 system and pH in the Chesapeake, and those that exist have focused on tributaries in the southern reaches of the estuary 27,28 . To address the coupling between acid-base chemistry and redox chemistry and its contribution to coastal OA, we sampled the water column repeatedly for several days within a deep basin of the main-stem bay in August 2013 and 2014, a time of peak hypoxia and anoxia, and these data were supplemented with an April 2015 (pre-hypoxia) study. In this paper, we report and explain the occurrence of a pH minimum at and above the oxic-anoxic boundary due to H 2 S oxidation. We further demonstrate how a combination of processes drives down pH and aragonite mineral saturation state, leading to CaCO 3 mineral dissolution in subsurface waters. Finally, we present a general geochemical model to explain why large eutrophic estuaries, exemplified by the Chesapeake Bay, are particularly vulnerable to the acidification stresses caused by the increase of anthropogenic atmospheric CO 2 and aquatic eutrophication and respiration.

Results
Acidification due to eutrophication-induced local CO 2 uptake. Partial pressure of carbon dioxide (pCO 2 ) in surface waters of the Chesapeake Bay exceeded 1600 μatm in its upper reach and was below the atmospheric level (~390 μatm) in the mid-bay for all three cruises between spring and summer ( Fig. 1a, b, Table 1, Supplementary Fig. 1, and Methods). The low pCO 2 was accompanied by high chlorophyll-a, a phytoplankton biomass proxy, and supersaturated dissolved O 2 for much of the year in the mid-and lower-bay , indicating net biological production fueled by high riverine nutrient loading. The low surface pCO 2 should lead to atmospheric CO 2 invasion  Table 1. The pCO 2 range was 340-590 µatm near CBL and 290-550 µatm near VIMS, with the lower ends represent incoming bay water during high tides and the high ends representing outgoing sub-estuarine waters during low tides and may contribute to water column CO 2 accumulation and acidification, particularly given atmospheric concentrations are ∼40% greater than the pre-industrial and bay has a long water residence time of 100 days 24 . While complete water column mixing and destratification occurs occasionally during storms 29,30 , smaller wind events more frequently mix water and chemical species down to middle depths (Fig. 2a, Supplementary  Figs. 6 and 7). Turbulence in the tidally driven bottom boundary layer will then mix the chemical species in the bottom water 31 .
We have estimated the air-to-water CO 2 flux and its impact on water column total dissolved inorganic carbon (DIC) and pH over the period of spring to summer (Table 1). To calculate the effect of DIC increase on bottom water pH decrease, we have modified the popularly used CO2SYS program to include H 2 S-HS − and NH 3 -NH 4 + species in the acid-base equilibrium calculations as the bottom water in August contains these reduced chemical species (Methods). The resulting bottom-water pH decreases (0.08-0.13 over the spring-summer period; Table 1) are significant when compared with pH decrease due to CO 2 uptake from the atmosphere in the open ocean (0.11). However the time scales of acidification due to local CO 2 uptake (months) are much shorter than open ocean uptake (decadal to centennial).
We also note that acidification induced by local CO 2 uptake is caused by both increased atmospheric CO 2 and coastal eutrophication. This is in sharp contrast with the CO 2 uptake in the open ocean where atmospheric forcing is comparable to the coastal ocean but biological CO 2 removal and physical mixing are less intense or frequent. Clearly, climate change, anthropogenic inputs, and natural processes have jointly altered the carbon cycle and stressed aquatic environments in the coastal zone. The spatial gradients of pCO 2 observed here and inferred net autotrophy are consistent with prior investigations using oxygenbased approaches to measuring primary production and respiration. Kemp et al. 25 concluded that the Chesapeake Bay was net autotrophic overall, but heterotrophic conditions (where respiration exceeded photosynthesis) prevailed in low-salinity regions where we measured supersaturated pCO 2 (Fig. 1a, [3][4][5]. Dissolved O 2 tends to be undersaturated in northern regions of the bay given high respiration rates associated with external loads of organic carbon 32 . These heterotrophic conditions gave way to a near balanced and an autotrophic metabolism in the mid-and lower-bay, leading to a mean, bay-wide net ecosystem production, which is consistent with strong under saturation of pCO 2 in these seaward regions (Fig. 1a, b and Supplementary Fig. 1). The mid-and lower-bay stations (CB3.3C and south) tended towards O 2 supersaturation during most months of the year, especially during the warmer months (Supplementary Figs. 2b, 4 and 5). Despite some interannual variability in the seasonal pattern of dissolved O 2 saturation, the years of 2013-2015 indicate similar seasonal patterns. Oxygen-based estimates of metabolism showed consistent surface-layer net O 2 production and bottom-layer net O 2 consumption, the rates of which were highly correlated 32 . Although oxygen-based methods could not be applied under oxygen-depleted conditions, independent measures of sulfate reduction (SR) in sediments, which dominated the benthic metabolism during warm months and led to significant sediment-water sulfide fluxes in the mid-bay 33,34 , clearly support the accumulation of sulfide observed in August 2013 (Fig. 2f).
Subsurface pH minimum due to oxidation of reduced chemicals. Repeated vertical profiles during both summers revealed a consistent pH minimum below the surface mixed layer at our focused study station, a deep site in the upper part of the mid-bay. Salinity profiles at this site ( Fig. 2a) combined with a time series of wind speed indicate a physical mixing event before our first sampling on August 9, 2013. Stratification quickly re-established when wind speed reduced and the wind direction switched from favoring mixing to favoring stratification ( Supplementary Fig. 6a, b). Total alkalinity (TA) and DIC were lower in the surface, but became higher in the bottom water (Fig. 2b, c). Dissolved O 2 was at saturation or supersaturation in the surface due to gas exchange and biological production and was not detectable below 10-15 m depths due to respiration (Fig. 2d).
On day 1, the mixed layer depth was still as deep as 15-18 m, but within 2 days, it shoaled to 10 m (Fig. 2a). Following this dynamic change, the O 2 penetration depth changed from about 15 m on day 1 to about 10 m on days 3-5 (Fig. 2d). Simultaneously, water column pH (25°C and NBS scale) decreased greatly over this period (Fig. 2e). For example, at the depth of 6 m, pH decreased from nearly 8.0 on day 1 to~7.5 on days 3-5. A pH minimum (7.35 ± 0.03) occurred at 11-13 m depth in the low O 2 zone (<10 μmol kg −1 ), below which pH increased slightly and then became constant at 7.45 ± 0.02. This pH minimum and the associated rapid pH decrease above it within a short period of <2 days have not been previously documented, although large pH changes were observed or expected in many strongly productive or stratified shallow water systems [35][36][37] . Such dramatic decreases in both O 2 and pH over just 2 days could likely put the biological system under stress 38 . This pH minimum was also observed in August 2014, though water column O 2 and pH profiles in August 2014 were relatively stable before and during the 5-days cruise, as wind speeds were less strong and less variable than those of the 2013 (Supplementary Fig. 7a, b), the oxycline (where O 2 decreases rapidly) and the pH minimum were even sharper and shallower, and bottomwater pH, DIC and TA were lower in 2014 than 2013 ( Supplementary Fig. 8).
We suggest the oxidation of reduced chemicals is responsible for the pH minimum in the low O 2 zone and the rapid pH decrease above it where declining O 2 gradients were steepest (Fig. 2e). The coupling between acid-base and redox chemistry is described by the following formula: Mn 2þ þ 0: Fe 2þ þ 0:25O 2 þ 1: In August 2013, because of the strong mixing event prior to our cruise, the total concentration of H 2 S was only 5 μmol kg −1 at the 20 m depth on day 1, but it rapidly increased to 30-40 μmol kg −1 on days 3-5 when the water column was restratified ( Fig. 2f) 39 . Oxidation of other reduced chemicals accumulating in the bottom water could also have contributed to the formation of the pH minimum. NH 4 + concentration measured near our site was 15-20 μmol kg −1 at 20 m ( Supplementary Fig. 9). Also, during days 3-5, bottom water [Mn 2+ ] and [Fe 2+ ] became as high as 7 and 2 μmol kg −1 , respectively 39 . When these reduced species (total concentration~60 μmol kg −1 ) were mixed upward into oxygenated water, they were oxidized, hydrogen ions were generated, TA was decreased and thus the water became more acidified (see Eqs. (1)-(4)). However, we recognize the oxidation of reduced species are often complex involving many intermediate steps and side products 40 and could have different H + production ratios.
It has been shown that oxidation of H 2 S by O 2 is sufficiently slow that H 2 S can be brought near to the surface during vigorous mixing events and lead to fish kills in coastal waters 20,36,41 . Similarly, ammonia oxidation is not instantaneous 12,17,41 . We suggest that the slow oxidation kinetics and rapid mixing facilitate the transport of reduced species and can subsequently result in acidification of the oxygenated near-and sub-surface waters, potentially resulting in a negative impact on aquatic organisms 38 . If, for example, one volume of bottom water of 60 μmol kg −1 of reduced chemicals is mixed with one volume of sufficiently oxygenated water, the resulting mixed water has the total concentration of the reduced chemicals halved to 30 μmol kg −1 , and eventually~60 μmol kg −1 of acid (or −ΔTA) would be generated due to the oxidation of the reduced chemical species (Eqs. (1)-(4)). Based on CO2SYS simulations, the predicted pH decrease due to these oxidation reactions can be up to 0.20 pH units in Chesapeake Bay waters, although other mixing ratios and incomplete reactions due to slow kinetics may generate less of a pH decrease (Methods). This pH decrease is substantial and is consistent with our observations ( Fig. 2e; also see a model simulation of TA, DIC, O 2, and pH evolving loci below). Note that while the size and location of this pH minimum may vary depending on the strength of the physical mixing and [H 2 S] in the bottom water, it occurs whenever bottom-water anoxia exists regardless of whether a prior severe mixing event has occurred as in our 2013 study, because moderate mixing occurs constantly in the bay (Supplementary Fig. 7).
In sediment porewater, a pH minimum was reported at and above the O 2 penetration depths as a result of oxidation of reduced chemicals, which diffused upward from deeper, anoxic depths 42,43 and was predicted by sediment diagenetic models 44,45 . Such a pH minimum was also seen in low O 2 waters of permanently stratified and anoxic deep basins including the Baltic Sea 13, 17 , the Black Sea 18,19 , the Framvaren Fjord 21 , the Hunnbunn Fjord 20 , and the Cariaco Basin 22 though no one has pointed out this phenomenon except Yao and Millero 21 who commented that "the low pH is difficult to explain". The pH minimum is an interesting feature that results from the decrease in TA:DIC ratio due to acid production during oxidation of reduced chemicals when encountering free O 2 due to vigorous physical mixing. To our knowledge, this is the first time that such a pH minimum has been reported and properly interpreted in the water column. We predict that the pH minimum should occur in all oceanic systems that have seasonally or permanently occurring oxic-anoxic boundaries, including the above mentioned cases as well as in the dead-end canals of Delaware Inland Bays 36 , Lake Grevelingen (the Netherlands) 12 , the Saanich Inlet 46 , and estuaries and bays elsewhere 17,20 . We further argue that the pH minimum is likely more dynamic in seasonally anoxic coastal systems than permanently anoxic deep basins, due to the shallower water depth and higher frequency of physical disturbances. Physical disturbances such as winds and tides occur regularly in the Chesapeake Bay ( Supplementary Fig. 7a, b) and other seasonally stratified coastal waters 12,36 . Therefore, their chemical and biological consequences, in the context of coastal OA and deoxygenation, deserve further attention.
Geochemical drivers and carbonate dissolution. To separate biological processes from physical mixing and to explore the biogeochemical control mechanisms in a broader context, we examine TA and DIC vs. salinity relationships at this site together with data from other areas of the bay and the river and offshore endmembers (Fig. 3a). Between the river and ocean endmembers, as expected, TA and DIC increased with salinity. However, at our focused study station, all subsurface and bottom-water samples were located well above the mixing lines, indicating net release of CO 2 and accumulation of DIC and TA. In addition, both DIC and TA data collected at the Chesapeake Biological Laboratory (CBL) dock, downstream of our focused study site at the lower end of the mid-bay, were also above the mixing lines. Those from the Virginia Institute of Marine Science (VIMS) dock, farther downstream near the bay mouth, however, showed the least enrichment relative to the conservative mixing lines. We also calculated the acid-base buffer factors from TA, DIC, and nutrients (PO 4 3− , H 2 S, and NH 4 + ) (Methods). It is clear the bay waters are poorly buffered as indicated by their much lower buffer factors compared to offshore waters here and elsewhere ( Fig. 3b; also see next section for definitions and explanations).
TA is usually a good conservative tracer of river-ocean mixing within an estuary because it is not influenced by CO 2 addition and removal. Because TA and DIC share a common major component (HCO 3 − ), deviations of DIC from the nearly conservative behavior defined by TA and salinity provide a measure of biological use or release of CO 2 47 . Bottom waters in the Chesapeake Bay, however, are conspicuously different from this general geochemical behavior normally encountered in oxygenated or moderately low oxygen environments 8,9,47 . DIC not only show a large enrichment against the conservative mixing line, but TA is also substantially enriched; with the excess DIC and TA reaching 275.3 ± 59.5 and 167.3 ± 54.2 μmol kg −1 respectively (Fig. 3a, Methods).
In any estuary, the most important internal sources of TA and DIC are aerobic respiration (AR), SR, and carbonate dissolution (CD) 12,13,21,45 (Fig. 4a, Table 2). Because each of these processes has a distinctly different ΔTA to ΔDIC ratio and involves a different pH change (Table 2), ΔTA:ΔDIC ratio and pH change become diagnostic of the geochemical processes. Based on the mixing line prediction, we can calculate the initial DIC and TA values at salinity (S) = 10 g kg −1 for surface water and S = 20 g kg −1 for bottom water. From the solubility constants we can also determine the initial concentrations of O 2 in S = 10 and 20 g kg −1 waters. Assuming the bottom water starts with a fully saturated dissolved O 2 , we can then derive DIC and TA generations and pH change for each step ( Table 2). In this poorly buffered water (Fig. 3b), the complete use of O 2 solely for AR would drive bottom-water TA and pH lower than the observed values (Fig. 4a, b). Sulfate reduction and CaCO 3 dissolution must then be invoked to explain the observed TA and pH. The effects of SO 4 2− reduction on TA and DIC can be estimated from the observed [H 2 S] ( Table 2) and the rest is made up by CaCO 3 dissolution (Methods). We envision that these three processes can occur either sequentially (Table 2) or simultaneously when O 2 , pH, carbonate mineral saturation state are sufficiently low. While the sequential pathway simulates the general patterns of the TA and DIC relationship (Fig. 4a) 6) and (7)). Most data were collected at the focused study site over a 5-day period, as marked in Fig. 1 (station 858). Those marked as TA other or DIC other were collected in waters between station 858 and the Chesapeake-Delaware Canal in the upper bay (Fig. 1). Data were also collected in the mid-and lower-bays at the CBL dock and the VIMS dock shortly after the main cruise  . 4b) reasonably well, it appears that CaCO 3 dissolution must have proceeded and occurred simultaneously with AR and SR as the simultaneous pathway simulates the observation better ( Fig. 4a, b) and as is justified by the very low aragonite carbonate saturation state once more than 50% of O 2 is consumed ( Table 2). From the above simulations, we conclude that up to~70% of the bottom-water TA production comes from CaCO 3 dissolution, which raises bottom-water pH from expected 7.25 to 7.45 and provides an important buffer mechanism in bottom waters (Methods). It has been reported that eutrophication led to lower pH in the polyhaline (S > 18 g kg −1 ) part of the bay between 1985 and 2008 26 . Based on our data and model simulations, we suggest that eutrophication in the bay has led to more O 2 consumption, SO 4 2− reduction, pH decrease, and dissolution of CaCO 3 shells and abiotic minerals in subsurface and bottom waters, consequently leading to possibly more TA and DIC export to the coastal ocean.
While shellfish calcification can represent a significant store of CaCO 3 in the Chesapeake Bay 26,48 , much of the CD present in the current study likely also comes from abiotic precipitation in surface waters, a mechanism noted before in the Loire estuary 49 . The extent of CD in deep waters estimated here could be supported by independently estimated CaCO 3 production in surface waters from deviations from conservative mixing ( Supplementary Fig. 10) 48 , which is consistent with the TA deficit observed here in low-salinity surface waters (Fig. 3a). While the precipitation may largely be driven by seasonal dynamics in primary production enhanced by estuarine eutrophication 49 , importantly, current and future increasing atmospheric CO 2 due to fossil fuel production may lower surface water carbonate saturation state enough to decrease mineral formation and thus delivery below the pycnocline. If so, the bay's deep water would have a reduced capacity to neutralize metabolically generated CO 2 , further enhancing eutrophication driven acidification.
Another important metabolic pathway is denitrification which uses NO 3 − as the oxidant for organic matter decomposition 50,51 . Note that [NO 3 − ] is generally low in the mid-bay (<1 μmol kg −1 ). However, denitrification is often coupled to nitrification at the sediment water interface. System-wide integrated denitrification rate has been estimated to be about 70 μmol m −2 h −1 in the Chesapeake Bay (summer time) 52 , although other estimations are lower. Taking this value as the upper end, we estimate that denitrification can contribute to a DIC production of up to 17 μmol kg −1 and TA production of up to 16 μmol kg −1 in a 10 m bottom-water column and over a 100-day period. This amount is only up to about 8% of the total TA production in the bottom water observed here. Finally, while organic matter decomposition using metal oxides as oxidants is important intermediate steps for biogeochemical cycles, the contributions to alkalinity production must be lower in the bay as recycled [Mn 2+ ] (<7 μmol kg −1 ) and [Fe 2+ ] (<2 μmol kg −1 ) 39 are much lower than the observed TA production in the bottom water; a conclusion similar to that derived in the Baltic Sea 13, 21 .

Discussion
The buffering capacity reflects the marine carbonate system's ability to resist changes in pH (or pCO 2 ) when DIC and/or TA are altered by physical and biogeochemical processes and when relevant thermodynamic constants are altered by temperature (T) and salinity (S) changes [53][54][55][56][57][58] . Mathematically, an aquatic system's ability to resist pH change can be deconstructed into its sensitivity to changes in T, S, DIC, and TA.
Here the first and second terms represent the effects of change in thermodynamic constants as a function of T and S. The third term reflects the pH change when DIC is added while keeping T, S, and TA constant and the fourth term reflects the pH change when a strong acid (H + or -ΔTA) is added while keeping T, S, and DIC constant. The slopes in the third and fourth terms are directly related to the buffer factors β DIC and β TA defined before 53, 54 with In estuarine conditions, because β DIC and β TA are similar in magnitude 54 (also see Fig. 3b), the overall contribution to acidification or pH decrease is largely decided by changes in DIC and TA during physical and biogeochemical processes (e.g., at constant T and S) and also by the initial buffering capacity (e.g., at variable T and S).
The pH and [O 2 ] relationship in Chesapeake Bay waters differs greatly from that observed in northern Gulf of Mexico (nGOM) waters (shaded line in Fig. 4b) 9 . It appears that Chesapeake Bay waters are more vulnerable to both anthropogenic CO 2 and biological induced acidifications because they have a lower buffering capacity than that of the offshore waters, in particular, in the nGOM as TA and DIC are lower in the Susquehanna River and US eastern margin waters than those of the Mississippi River and nGOM seawater (Figs. 3b, 4b and 5a) 9, 59, 60 . However, our simulations and those of the previous studies 12, [53][54][55][56][57] suggest that lower buffering capacity itself does not necessarily lead to low pH (Fig. 5a); rather, it allows a much greater pH decrease when other sources of CO 2 or strong acids are added (Fig. 5b). Similar amounts of AR (Table 2) would lead to a pH decrease of only 0.4 units in the strongly buffered nGOM waters whereas a larger decrease of nearly 0.8 units would occur in the poorly buffered Chesapeake Bay waters at the present day conditions (at S = 34 and 20 g kg −1 , respectively, Figs. 4b and 5c).
While the OA signal due to CO 2 uptake in the open ocean regions is similar across middle and lower latitudes, the manifestation of this anthropogenic CO 2 signal through ocean-river mixing in estuaries is dependent on the river TA and DIC values, which are highly variable among the world's rivers 61 , and whether CO 2 is also introduced via microbial respiration 11,60 . Due to the very high river TA and DIC and the resulting strong buffering capacity over the entire salinity range in the Mississippi River impacted coastal waters, pH change due to OA is proportional to the open ocean OA source signal and salinity and decreases toward zero salinity (Fig. 5a) 9,60,62 . In the Chesapeake Bay where average river TA, DIC, and buffering capacity are low, however, the oceanic OA signal is amplified in the low and middle salinity zone. Here the combination of reduced buffering capacity (with decreasing salinity) and a still sufficiently strong open ocean OA signal generates a minimum buffer zone 60 and thus a Maximum Estuarine Acidification Zone (MEAZ) (Fig. 5c). The existence of a MEAZ and its salinity range depend not only on the river TA value, but also the TA:DIC ratio 60 . When CO 2 addition from AR increases from 0 to 100 and finally to~200 μmol kg −1 (or 0 to roughly half or to a full O 2 consumption depending on the salinity and temperature), the minimum buffer zone shifted from salinity~4 to~13 (Fig. 5a) and finally to~23 (Fig. 5b) 60 . Note that local CO 2 uptake, carbonate mineral dissolution and SR are not included in this discussion (Fig. 5) and would further modify the estuarine buffering capacity as they would modify the TA to DIC ratio in estuarine waters ( Table 2, Figs. 3b, 4a, b).
Below we further discuss the effects of anthropogenic CO 2 and biological CO 2 and acid additions on estuarine pH buffering capacity 53,54 . The marine carbonate system has a minimum buffering or maximum pH change point when DIC increases approximately equal to that of TA (or TA:DIC ≈ 1) where [CO 2 ] = [CO 3 2− ] + [B(OH) 4 − ] (if we ignore all other weak acid-base species). At this point, any addition or removal of CO 2 or acids will result in a maximum pH decrease or increase. Because DIC is slightly higher than or nearly equal to TA in rivers 47,62,63 and is lower than TA in seawater, there may exist a small crossover of DIC and TA at the very low-salinity zone. A peculiar pH minimum occurring in the low-salinity zone of estuaries is related to this mixing feature as was noticed a long time ago in both closed 64 and open 65 system simulations. Furthermore, how CO 2 is added to the estuarine waters affect how the crossover point will move. First we contend anthropogenic CO 2 does not directly add to the high pCO 2 river water but is mixed into the estuary via river-ocean mixing 60 . In contrast, respiratory CO 2 is nearly equally added to the bottom water based on O 2 consumption regardless of the mixing index or salinity (except that O 2 solubility increases when salinity decreases, but it is a small correction). In Fig. 6, we summarize several scenarios illustrating how the crossover point of the TA and DIC to salinity lines or the point of TA:DIC ratio = 1 moves along the TA-salinity line. Adding anthropogenic CO 2 to the seawater endmember would move this crossover point to only a slightly higher salinity. Adding biological CO 2 (for example 100 μmol kg −1 ) to both the river and ocean endmembers would, however, shift the DIC line to a much higher position (parallel to the original line) creating a crossover point located at a salinity substantially higher than the original one. Finally the combined effect of anthropogenic CO 2 and biological CO 2 from respiration moves the crossover point to an even higher salinity. These crossover points are consistent with the progressive shift of the minimum buffer factor (β DIC ), the pH minimum, and the maximum acidification zone (-ΔpH) presented in Fig. 5. However there appears a difference in the location (salinity) between the minimum buffer factor (β DIC ) and the maximum acidification (-ΔpH). This is because ΔpH represents the accumulative pH change between an end point and a beginning point while the buffer factor reflects the pH sensitivity at a specific point if additional DIC is added to the system. In summary, large estuarine water bodies, exemplified by the Chesapeake Bay, are particularly vulnerable to the anthropogenic CO 2 and eutrophication-induced coastal OA. In this paper we emphasize that subsurface H 2 S oxidation (~0.2 pH units) and local surface CO 2 uptake (~0.1 pH units) work together with known organic matter respiration and the open-ocean OA source signal to drive substantial acidification and CaCO 3 dissolution in estuarine subsurface waters. Currently, acidification due to CO 2 input from AR (up to 0.8 pH units) in the bay exceeds that from the atmospheric CO 2 increase in the open ocean (~0.1 pH units in surface waters and the signal is mixed into estuaries proportional to salinity) and local estuarine uptake; but towards the end of this century the latter will approach or exceed the former and the synergy between them will also increase. In addition, future increasing atmospheric CO 2 due to fossil fuel production may lower carbonate saturation state enough to decrease mineral formation in surface water and thus delivery below the pycnocline, where we have currently demonstrated that CaCO 3 dissolution offsets a significant proportion of the metabolic CO 2 effect on acidification. We further recognize that natural and anthropogenic acidification mechanisms most relevant to estuarine acidification are characterized by various time scales. They range from nearly instantaneous for acid-base equilibrium 66 , to minutes for CO 2 hydration 66 , and to minutes to hours for H 2 S oxidation 36,41,67 . In contrast, time scales for physical mixing are on the order of tidal or less, to daily and seasonal 30, 31, 68 while local CO 2 uptake from the atmosphere and its accumulation in the water column and acidification of the bottom water as well as pelagic and benthic respirations operate over tidal to seasonal scales 12,13,35,37,69 . Although anthropogenic changes in external forcing due to variability in river and ocean endmembers may also have a seasonal component, decadal and centennial variation is more important 5,35,70,71 . This mosaic of processes with different time and space scales poses a great challenge in our ability to understand and predict coastal OA.

Methods
Site and cruise descriptions. The Chesapeake Bay is the largest estuary in the US. The August 2013 survey started from the upper estuary near the Susquehanna River mouth (Fig. 1a). The upper and middle Chesapeake Bay were surveyed during 9-14 August 2013 by RV Hugh R. Sharp. The water column survey focused mainly at one site south of the Bay Bridge (38°58.8 N, 76°22 W), where a field study  a difference between b and a). The red arrows represents the Maximum Estuarine Acidification Zone (MEAZ). Note that β DIC presented here does not match that from the real data (Fig. 3b) at low salinities as a stratified and O 2 depleted state does not exist for the very low salinity part of the Chesapeake Bay  Fig. 6 Estuarine TA or DIC mixing lines between the river and ocean endmembers with additional CO 2 added in from anthropogenic sources and respiration. Line 1 is river-ocean mixing of TA, line 2 is river-ocean mixing of DIC, line 3 is mixing DIC plus anthropogenic CO 2 addition, line 4 is mixing DIC plus respiratory CO 2 addition, and line 5 is mixing DIC plus both anthropogenic and biological CO 2 additions. The crossover points are river-ocean mixing lines of TA with A river-ocean mixing lines of DIC, B mixing plus anthropogenic CO 2 , C mixing plus CO 2 from respiration, and D mixing plus CO 2 from both B and C of redox chemistry 39 and water column inorganic carbon and pH were carried out. We consider this site as the up end of the mid-bay. We repeatedly sampled the water column at high slack tide and low slack tide. During this 5-day survey, an excursion was made south to the middle bay near Solomons Island. After the completion of the cruise, we set up two 24 Fig. 1). We also conducted a spring survey (11-16 April 2015) to get an initial condition before the hypoxia season (Fig. 1b).
Sample and analytical methods. Salinity, temperature, and O 2 were obtained from the CTD Rosette system. Total sulfide (H 2 S and HS − ) was determined by voltammetry using solid state Au/Hg electrodes 36,39,43 . Surface water partial pressure of CO 2 (pCO 2 ), position, temperature, and salinity information were measured underway while the ship was sailing or anchored by pumping surface water from under the ship to the shipboard laboratory using an underway pCO 2 system 72 . TA and DIC water samples were taken from Niskin bottles and were preserved and stored in 250 ml borosilicate glass bottles with 100 μl saturated HgCl 2 solution 73 . TA and DIC samples were stored at refrigerated temperature (~5°C ) before being measured (within 4 weeks). TA samples were measured by opencell Gran titration with a precision better than ±0.1% using an Apollo Scitech Seawater Total Alkalinity titration system 73 . DIC samples were analyzed by adding phosphoric acid into sample waters to release CO 2 , which was measured by an infrared CO 2 analyzer (LI-COR 7000) with an overall precision of ±0.1% using an Apollo Scitech DIC Analyzer 73 . Both TA and DIC measurements were quality controlled by Certified Reference Materials from Andrew Dickson of the University of California at San Diego. pH samples were taken by the same Niskin bottle and were measured by an Orion Ross glass electrode within 1 h after the water temperature was stable in a 25.0 ± 0.1°C thermal bath on the research vessel. The electrode was calibrated against three NBS (NIST) standards and pH values are reported in NBS scale and at 25°C. Note that pH values in NBS scale are about 0.1 pH unit higher than those reported in total H + scale (pH T ) elsewhere.
Uncertainty in determining DIC and TA enrichment. We averaged all subsurface and bottom waters (S > 11 g kg −1 ) to derive a DIC enrichment of 275.3 ± 59.5 and a TA enrichment of 127.3 ± 54.2 μmol kg −1 , with respect to their expected conservative behaviors. The TA enrichment is underestimated because of a technical challenge caused by HgS precipitation, which releases H + when HgCl 2 was used to stop microbial activity 18 , and/or oxidation of H 2 S and NH 4 + during sample storage and/or analysis (which also generate H + ). TA values calculated from DIC and pH analyzed onboard (neither is subject to the sample preservation and storage problems) agree well with measured TA except in the bottom waters ( Supplementary  Fig. 11a), and the disagreement increases as [H 2 S] increases ( Supplementary  Fig. 11b). The internal consistency analysis suggests that TA reduction due to sample preservation and storage is 40 ± 20 μmol kg −1 ; thus, the most likely TA enrichment in the bottom water is (127. Calculation of air-sea CO 2 flux and the DIC increase. Water surface pCO 2 was measured every 1-1.5 min with calibrations every 6-12 h (August 2013, August 2014, and April 2015). Atmospheric pCO 2 values were also measured every 2-4 h during these cruises. Both atmosphere and water CO 2 values were measured in a dry condition (xCO 2 ) and were converted to pCO 2 in 100% water saturated conditions inside the equilibrator (pCO 2(eq) ) by considering water vapor pressure: where P b is barometric pressure and P weq is water vapor pressure in the equilibrator. For water data, the pCO 2(eq) is further converted to estuarine surface water pCO 2 (pCO 2(water) ) by considering temperature changes between the surface water and the equilibrator through the following equation 72 : where SST is sea surface temperature (°C) and T eq is temperature in the equilibrator. Our measured atmospheric xCO 2 values were also converted from dry condition to near sea surface wet condition (pCO 2(air) ) by Eq. (10): Here P b is barometric pressure and P w is water vapor pressure at the sea surface. Each pCO 2(water) and its corresponding pCO 2(air) were used to calculate the gas exchange flux (F CO2 ) between atmosphere and water by Eq. (11).
where k represents the gas transfer velocity and K 0 is the solubility of CO 2 74 . We adopted Ho et al. 75 as the gas transfer velocity and an ensemble of gas transfer parameters to evaluate the uncertainty range following the previous practice 72,76 . Finally, the coefficient C 2 corrects the non-symmetrical distribution of wind 76 . A negative air-sea CO 2 flux means an uptake of atmospheric CO 2 for the water.
Over at least a 100-day water residence period (from May to August) and over a water column of 20 m, this CO 2 flux can be converted into an increase in DIC of, 4.3-6.7 (mmol m −2 d −1 ) × 100d/20 m ≈ 21.5-33.5 mmol m −3 or 21.2-33.1 μmol kg −1 (here a density of 1012.09 kg m −3 is used). We used the entire water column rather than the surface mixed layer because the main concern here is how local CO 2 uptake, via internal mixing, contributes to acidification of the especially vulnerable bottom waters.
Calculation of pH decrease due to local CO 2 uptake. With H 2 S included, the calculation of pH (in NBS scale and at 25°C) decrease was performed using the modified CO2SYS program. Note another program, AquaEnv, also has such a capacity 77  Then, we subtracted the summer DIC by 21.2-33.1 μmol kg −1 (=1912. 6-1900.7 μmol kg −1 ) to calculate a new pH (7.557-7.604). Thus, the pH decrease by an increase of DIC derived from local uptake of atmospheric CO 2 is 0.081-0.128 pH unit over the entire period from spring to summer.
pH decrease due to oxidation of reduced chemicals. We used day 4 data with S = 15.145 g kg −1 , T = 25.36°C, depth = 12.54 m, DIC = 1767.7 μmol kg −1 , [H 2 S] = 2 μmol kg −1 , and pH = 7.354 to calculate a TA = 1698.8 μmol kg −1 . Then, we subtracted a 30-60 μmol kg −1 from TA to calculate a new pH (7.152-7.246). Thus, the pH decrease by a 30-60 μmol kg −1 of TA reduction is 0.108-0.202 pH units (represented by the purple arrows in Fig. 4b). The modified version of CO2SYS was used for all the CO 2 and pH calculations. Note, adding a <3 μmol kg −1 of T-PO 4 would only lead to <0.005 pH unit decrease in the calculation. Thus, its infleunce is ignored here.
Modification of the CO2SYS program. The modifications were done on the Excel version 2.1 of the program 78 , which is available for download from CDIAC (http:// cdiac.ornl.gov/ftp/co2sys/). A Matlab version is available from the corresponding author. In addition to the total Phosphate and Silicate, the program now accepts the total NH 3 and total H 2 S in µmol (kg of SW) −1 . The contribution of each to the alkalinity is given by: and where K NH3 and K H2S are the dissociation constants of ammonium (NH 4 + ) and hydrogen sulfide (H 2 S).
The dissociation constant for NH 4 + was taken from Clegg and Whitfield 79 and is valid for S = 0-40 g kg −1 and t = −2 to 40°C (note ref. 21 essentailly provided the same constant). The constant for H 2 S was taken from Millero et al. 80 and is valid for S = 0-40 g kg −1 and t = 0-35°C. When the pressure is not zero, a correction is applied according to Millero 81 . A comparison with AquaEnv under [H 2 S] <50 (or 300) μmol kg −1 shows a good agreement of calculated pH (in free scale) within 0.0003 (or 0.0026) from known TA and DIC. We have further tested the calculations with waters containing high concentrations of H 2 S and posted this modified version of the CO2SYS program on the CDIAC website for public access 82 . total of 195.9 μmol kg −1 (i.e., observed 167.3 μmol kg −1 plus expected −28.5 μmol kg −1 ) beyond conservative mixing. We estimate TA increase from SO 4 2− reduction as 80.6 μmol kg −1 from the total concentration of H 2 S (35 μmol kg −1 ) 39 by the following equation: where 2 × [H 2 S] represents an equal amount of HCO 3 − and HS − production during SO 4 2− reduction and 16/53 × [H 2 S] represents NH 3 production (and contribution to TA) based on stoichiometry (see Table 2). Then the TA generated from CaCO 3 dissolution must be as high as 115.3 ± 20.0 μmol kg −1 by the difference (195.9-80.6) and contributes up to 70% of total amount of TA production. The amount of DIC production following these steps is 305.6 ± 10.0 μmol kg −1 . This is within the uncertainty of the observed value of 275.3 ± 59.5 (Fig. 3). The 10% difference (30 μmol kg −1 ) can be explained either by TA increase due to organic matter respiration using nitrate (denitrification) and metal oxides 13,45 and/or the deviation of C/N ratio from the Redfield ratio 13,84,85 as well as probably organic alkalinity contribution 86 . Indeed if a lower C to −O 2 ratio (106/154 = 0.688) given in ref. 84 is used, the produced DIC would be close to the observation.
From the resulting pH and Ω arag ( Table 2), it is clear when O 2 is partially consumed aragonite mineral becomes undersaturated (starting at 75% O 2 saturation) and CaCO 3 dissolution can proceed together with AR. To simulate the observed DIC (measured) and TA (calculated from DIC and pH) data, we assume the dissolution does not occur until a sufficiently low Ω arag in waters and that the first 20% of CaCO 3 dissolution occurs before or at DO = 33%. Then the second, third and fourth 20% of the CaCO 3 dissolution occurs before or at DO = 16%, 8%, and 0%, respectively. The last 20% of the CaCO 3 dissolution occurs together with SR.
Simulation of pH changes. We present here the calculation methods for Fig. 5a, b. Although each term in Eq. (5) may be derived analytically, in this paper, we obtain the overall pH change, ΔpH (presented in Fig. 5c), numerically using the updated CO2SYS program. Conditions used are given below.
We assume a present day atmospheric dry CO 2 fraction (xCO 2 ) as of 396.2 ppm and the corresponding water pCO 2 = 384.0 μatm at 25°C and salinity = 36.0 g kg −1 . xCO 2 is set to 281.3 ppm for the pre-industrial era and 798.0 ppm for year 2100.
For the Chesapeake Bay simulation, we take TA from the offshore water at 74. Buffer factor calculation. We calculate these buffer factors (Fig. 5) following the analytical formula provided by Egleston et al. 54 with a typo corrected 87 . Specifically, we extract out species concentrations and thermodynamic constants from the updated version of CO2SYS. We also compare results with those calculated with two other programs AquaEnv and SEACARB (http://CRAN.R-project.org/ package=seacarb) and the agreement is reasonable (within 8 μmol kg −1 or 3%) as is expected. Although H 2 S and NH 4 are not included in the analytical equations, as our equilibrium calculation already include these acid-base species and [H 2 S] and [NH 4 ] are not high in the bay, the calculated results are similar to those from the AquaEnv which includes fully these reduced chemical species. Finally, in Fig. 3b, for the real system buffer factors, we directly use AquaEnv. Given the low [H 2 S] in the bay, with or without H 2 S have only resulted in a minor difference in buffer factor calculation.
Computer program availability. The modified CO2SYS program on Excel version 2.1 is available for download from CDIAC (http://cdiac.ornl.gov/ftp/co2sys/). The Matlab version is available from the corresponding author upon reasonable request.
Data availability. All data are available from the corresponding author upon reasonable request and will be deposited at the US National Centers for Environmental Information (https://www.nodc.noaa.gov/oceanacidification/).