Influence of manganese cycling on alkalinity in the redox stratified water column of Chesapeake Bay

The alkalinity dynamics in coastal environments play a crucial role in controlling the global burial of carbonate minerals and the ocean's capacity to sequester anthropogenic CO₂. This study presents results from high vertical resolution profiles obtained during two summers in the temperate Chesapeake Bay estuary, enabling detailed investigation of carbonate dynamics over salinity and redox gradients, along with measurement of the speciation of most redox-sensitive elements. Under oxygen-rich conditions, carbonate dissolution, primary production and aerobic respiration explain the evolution of total alkalinity (TA) versus dissolved inorganic carbon (DIC), once adjusted for fresh and oceanic water mixing. A significant flooding event in 2018 promoted carbonate dissolution. In oxygen-depleted waters, we observed a previously unreported 2.4 mol increase in DIC per 1 mol of TA production, which was consistent over the 2 years. Stoichiometric changes suggest that MnO₂ reduction followed by Mn carbonate precipitation is responsible for this characteristic carbonate signature, likely produced in sediment pore water and then transferred to the water column along with other by-products of anoxic respiration at the onset of summer. Our findings highlight the critical role of Mn in alkalinity dynamics in the Chesapeake Bay and potentially other river-dominated environments where it can limit H₂S oxidation to ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ and promote sulfur burial.

1 Introduction

About 30 % of anthropogenic CO₂ emissions are rapidly trapped by dissolution in the ocean as dissolved inorganic carbon (DIC), which is dominated by bicarbonate ions (${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$; Friedlingstein et al., 2019). At the century timescale, atmosphere–ocean exchanges result in oceanic ${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$ enrichment not associated with a cationic enrichment, as measured by total alkalinity (TA). This is in contrast to silicate or carbonate weathering, which preponderates at timescales of thousands to millions of years and increases TA simultaneously with DIC (Urey, 1952). This cationic deficiency corresponds with an excess of proton release that can be tracked by an indicator such as TA minus DIC (Xue and Cai, 2020) and informs about the reduction of the CO₂ buffering capacity of seawater (Gustafsson et al., 2019). The proton excess is only balanced in the deep ocean by the production of Ca²⁺ from carbonate dissolution, a process named chemical carbonate compensation (Boudreau et al., 2018). However, in shallow waters, that accounts for two-thirds of global buried carbonate (Smith and Mackenzie, 2016); carbonate precipitation largely predominates over dissolution, and other localised processes may constrain carbonate dynamics (Borges et al., 2006; Lohrenz et al., 2010). For example, some calcifying species may slow down their carbonate precipitation in case of pH decrease (so-called biological carbonate compensation; Boudreau et al., 2018). As another example, in estuaries, the seasonality of river flow, temperature and continental erosion modulates CaCO₃ dissolution (e.g. Su et al., 2020b), atmospheric CO₂ exchanges (e.g. Borges et al., 2018), respiratory activity (e.g. Abril et al., 2003) and transfer of carbonate particles (Meybeck, 1987; Middelburg et al., 2020).

The global trends of human migration towards littoral areas and global warming favour eutrophication and a decrease in oxygen levels in coastal water (Breitburg et al., 2018; Rabalais et al., 2014). Sometimes driven by recurrent natural process (e.g. Gupta et al., 2021), exceptional events (e.g. Hulot et al., 2023), stratification reinforced by global warming (e.g. Meire et al., 2013) or anthropogenic nutrient loading (e.g. Hagy et al., 2004; Carstensen et al., 2018), coastal deoxygenation enhanced the possible build-up of anoxic conditions to build up as a permanent or seasonal feature. In anoxic environments, the drivers of the carbonate cycle change as the terminal electron acceptor for organic matter remineralisation shifts from O₂ to Mn oxides, Fe oxides and then sulfate. Notably, the consumption of sulfate, followed by sulfur burial, produces an significant alkalinity flux toward coastal waters, accounting for approximately 10 % of the global alkalinity input to the ocean (Middelburg et al., 2020). Fe and Mn oxides play a critical role in sulfur burial, as they prevent ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ regeneration and lead to the production of S⁰ or FeS instead (Findlay et al., 2014; Avetisyan et al., 2021), leaving a specific fingerprint in the TA and DIC effluxes.

While TA controls the CO₂ buffering capacity of the ocean, riverine input of carbonate to the ocean is poorly constrained (Middelburg et al., 2020), and only rare publications take into account the estuarine transformations of the carbonate species (e.g. Su et al., 2020a; Abril et al., 2003) furthermore in the context of oxygen depletion (e.g. Abril et al., 2004). To better constrain the carbonate cycle in oxygen-depleted estuaries, we sampled a stratified water column in the Chesapeake Bay 11 times over two campaigns with a high vertical resolution (down to 10 cm). This protocol allows a precise description of carbonate dynamics over a redox gradient along with the measurement of the speciation of most redox-sensitive elements. Such sampling illustrates carbonate chemistry on transitioning from oxygenated waters to waters devoid of oxygen, as usually only encountered in sediments or in anoxic lakes or seas (e.g. the Black Sea). The original observed changes of alkalinity versus dissolved inorganic carbon changes are interpreted based on typical geochemical reactions occurring along the redox gradient.

2 Materials and methods

2.1 Sampling

During two sampling campaigns from 3 to 9 August 2017 and 28 July to 3 August 2018, 11 profile casts were conducted in a single station in the Chesapeake Bay with a water depth of 25 m (station 858; $\mathrm{38}{}^{\circ }{\mathrm{58.54}}^{\prime }$ N, $\mathrm{076}{}^{\circ }{\mathrm{22.22}}^{\prime }$ W). The Susquehanna River is the main tributary of the bay representing on average two-thirds of the fresh water input (Zhang et al., 2015). Despite the similar season, the two campaigns occurred at very different river flow, with about 850 m³ s⁻¹ in 2017 versus 8500 m³ s⁻¹ in 2018 due to release of flood waters from the Conowingo Dam. The August 2018 condition corresponds to flooding which occurs on average every 3.5 years (return period of 3.5 years, USGS survey).

Each CTD (conductivity, temperature, depth) cast was performed during low or high tide slack periods. An oxygen sensor (Clark electrode, SBE Inc.; detection limit of 1 µM) and fluorescence sensor (ECO-FL Fluorometer, WET Labs) were part of the CTD Rosette to take measurements during sampling. Also, a submersible all-plastic pump profiler was attached with the pump near the sensor orifices, allowing measurement and sampling at a resolution of a few centimetres over 25 m water depth. Water was pumped to the deck within 1 min, and water passed through a flow through voltammetry system measuring O₂, Mn(II), Fe(II), organically complexed Fe(III), FeS clusters, H₂S and polysulfides continuously (Hudson et al., 2019). When redox interfaces were identified, samples were filtered through an acetate cartridge filter (pore size 0.45 µm) for pH and inorganic carbon parameters, which were processed onboard within a few hours after sampling in order to conserve chemical speciation. The pump profiler system was cleaned with deionised water (18 MΩ) onboard the deck of the ship after deployment. No coating effects were observed with the pump system.

2.2 Discrete measurements

For each sample, all redox species were determined in the through flow voltammetry system using cyclic voltammetry with a 100 µm diameter Au/Hg amalgam PEEK microelectrode prepared according to Luther et al. (2008) connected to a DLK-60 electrochemical analyser from Analytical Instrument Systems Inc. The detection limit of this method is 0.2 µM for sulfide and polysulfides. Discrete samples for the determination of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, Fe and Mn species were filtered through nylon luer-lock syringe filters (Millipore, 0.20 µm). Iron was measured on both bulk and filtered samples using the ferrozine method (Stookey, 1970): after HCl acidification (for Fe(II)) and an optional reduction step (for Fe(III)+Fe(II)) with hydroxylamine hydrochloric (final concentration 0.7 M for 1 h), ammonium acetate (final concentration 0.5 M) and ferrozine (final concentration 1 mM). Absorbance at 562 nm was read with a diode array spectrophotometer (Hewlett Packard 8452B). Limit of detection is 100 nM for Fe(II) and Fe(III) with a 1 cm cell. Shipboard nitrite determination was performed using the method of Grasshoff (1983). To 25 mL of sample, 0.5 mL of 58 mM sulfanilamide in 10 % $v/v$ HCl and 0.5 mL of a 4 mM N(1-naphthyl)ethylene diamine hydrochloride solution were added. Samples with added reagents were shaken and left to sit for 15 min, followed by UV–Vis analysis at 540 nm using a 10 cm cell to increase detection limits. Calibration curves were constructed using sodium nitrite. Limit of detection is 10 nM for ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$.

Dissolved manganese was determined by displacement of a Cd(II)-porphyrin complex with Mn(II) to form the Mn(II)-porphyrin complex (Ishii et al., 1982). Mn(III) species were identified based on slower reactivity with the Cd complex (Madison et al., 2011) as modified in Thibault de Chanvalon and Luther (2019). Alternatively, Mn(III) species were identified after HCl treatment (down to pH=1.5) followed by filtration in order to flocculate and eliminate the dissolved manganese bound to humic material by filtration (Oldham et al., 2017b). Limit of detection is 50 nM for Mn(II) and Mn(III) in a 1 cm cell. MnOx was measured on 20 mL samples of suspended material retained on 0.2 µm filters by the Leucoberbelin blue (LBB) method (Jones et al., 2019). Four millilitres of a reagent solution ([LBB]=78 µM, [acetic acid]=14 mM) react with the filter, and the absorbance is read at 624 nm. KMnO₄ was used to calibrate the LBB method, which allows the calculation of the electron equivalents obtained from particulate MnOx. Results are given as MnO₂ equivalent, with a limit of detection of 0.1 µM and an uncertainty below 5 %.

The DIC samples were preserved in 250 mL borosilicate glass bottles with 50 µL saturated HgCl₂ solution. The total alkalinity (TA) samples were not poisoned to prevent HgS precipitation and H⁺ release in anoxic and low-salinity waters (Cai et al., 2017). Then TA was analysed by Gran titration in an open-cell setting (AS-ALK2, Apollo SciTech) within 24 h of collection (Cai et al., 2010). The DIC samples were measured by a nondispersive infrared analyser (AS-C3, Apollo SciTech) within a week (Huang et al., 2012). The precision for DIC and TA was about 0.1 %. Both DIC and TA measurements were calibrated against certified reference materials (CRMs Batch 163 and Batch 173 provided by Andrew Dickson of the Scripps Institution of Oceanography). The pH samples were measured onboard at 25 ^∘C within 1 h of collection using an Orion ROSS glass electrode and calibrated with NIST standard buffers. The pCO₂, calcite saturation and TA were calculated from measured DIC and pH via the CO2SYS program using Cai and Wang (1998) constants. The measured TA was found to be highly correlated to the calculated TA (r²=0.995 and 0.998, slope=0.995 and 1.017 for the 2017 and 2018 campaign respectively), and their difference was always below 45 µM with an average of 7.5 µM for 2017 and of 22.2 µM for 2018. These results suggest a low contribution of non-carbonate species (e.g. nitrite, ammonium or organic matter; Cotovicz et al., 2016), and measured TA was used for the interpretation.

2.3 Models of biogeochemical process on TA and DIC

2.3.1 Identification of biogeochemical processes from scatter plots

This section identifies the required conditions to interpret a scatter plot of two species in terms of chemical reaction stoichiometry: the “reaction-driven” approximation. In case of turbulent diffusion mixing (sometimes called eddy diffusion) in only one direction (no lateral input), at steady state, on a portion of space where occurs one chemical reaction, the changes of concentration, C and D, of two species can be described by Eq. (1) (see Sect. A1 for more details):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}x}}={\displaystyle \frac{{\mathit{\alpha }}_C}{{\mathit{\alpha }}_D}}{\displaystyle \frac{\mathrm{d}D}{\mathrm{d}x}}+G,\end{array}$$

with α_C and α_D the respective stoichiometric coefficient, x the direction in the space and G an unknown constant. Here, C and D are concentrations of any elements produced or consumed during the reaction and can be TA, DIC or a linear combination with salinity. When G is null, Eq. (2) is valid and the data points would be distributed along a straight line in a C versus D scatter plot.

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}D}}={\displaystyle \frac{{\mathit{\alpha }}_C}{{\mathit{\alpha }}_D}}\end{array}$$

Theoretically, Eq. (2) is valid in at least three cases. First, it is valid if the contribution from one endmember is negligible compared to the timescale of reactions, which corresponds also to a system with only one endmember. Second, it is valid when the slope of $\mathrm{\Delta }C/\mathrm{\Delta }D$ between the endmembers equals ${\mathit{\alpha }}_C/{\mathit{\alpha }}_D$, which is achieved, in particular, if one endmember was previously generated from the other endmember by a chemical reaction with similar stoichiometry. And third, it is valid when the two endmembers have similar concentrations (see Sect. A2), as in this study where the excess of TA (TAex) and the excess of DIC (DICex) are calculated by linear combination with salinity to be equal to zero for the upstream and downstream endmembers. Thus, if only one reaction occurs in between the endmembers, the slope $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ would correspond to the stoichiometry of ${\mathit{\alpha }}_{\text{TA}}/{\mathit{\alpha }}_{\text{DIC}}$. The oceanic endmember for both August campaigns was the one proposed by Su et al. (2020a). The oceanic endmember varies mainly with the season (Cai et al., 2020), and a maximal change of 50 µM results in 5 % uncertainty on the slope of the mixing line. For the upstream estuary, larger variations exist mainly due to changes of river discharge and water temperature (Meybeck, 2003; Joesoef et al., 2017), and a one-off endmember has to be determined by fit with the in situ measurements at the lowest measured salinity. Between the two campaigns, the upstream endmember changed by 77 µM, generating 5 % of change on the slope (see Fig. A1). The upstream endmember is not a river endmember (Su et al., 2020a) but corresponds to a salinity above 1.5, avoiding any interpretation for biological activity in the fresh water part of the estuary (Meybeck et al., 1988). However, it corresponds to a larger water mass pool, less sensitive to short term changes, with a residence time being higher than 240 d in the Chesapeake upstream part (Du and Shen, 2016), and thus is more likely to satisfy the condition of steady state.

The relative uncertainty on $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ is equal to the sum of the relative uncertainty of ΔTAex and ΔDICex. Posing ΔTA as the change of TA measured, ΔS as the change of salinity and sml_TA as the slope of the mixing line for TA, we have $\mathrm{\Delta }\text{TAex}=\mathrm{\Delta }\text{TA}-\text{sml\_TA}\times \mathrm{\Delta }S$. Uncertainty on ΔTA and ΔS are negligible to the relative uncertainty of the slope of the mixing line. Posing δ(x) as the uncertainty on x, we get Eq. (3) that describes the fact that the uncertainty is much lower on ΔDICex than on DICex because most of the error associated with the calculation of the endmember is cancelled when calculating the difference of DICex on two points with close salinity:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\mathit{\delta }(\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex})}{\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}}}& ={\displaystyle \frac{\mathit{\delta }\left(\mathrm{\Delta }\text{TAex}\right)}{\mathrm{\Delta }\text{TAex}}}+{\displaystyle \frac{\mathit{\delta }\left(\mathrm{\Delta }\text{DICex}\right)}{\mathrm{\Delta }\text{DICex}}}\\ & ={\displaystyle \frac{\mathit{\delta }\left(\text{sml\_TA}\right)}{\text{sml\_TA}}}+{\displaystyle \frac{\mathit{\delta }\left(\text{sml\_DIC}\right)}{\text{sml\_DIC}}}\\ & =\mathrm{0.1}.\end{array}\end{array}$$

However, in a stratified water column, not just one but several successive reactions occur, limiting the validity of Eq. (1) to each reactional stratum. The general case is not straightforward to solve, but in the particular case where the C versus D plot represents a straight line between two endmembers with different concentrations, the previous discussion of Eq. (2) indicates (second case) that one endmember would have been previously generated from the second by a chemical reaction with similar stoichiometry. Thus, the depths corresponding to the straight line define a reactional stratum characterised by a constant ${\mathit{\alpha }}_C/{\mathit{\alpha }}_D$ and delimitated by two local endmembers maintained in steady state by chemical reactions with similar stoichiometry than the one that produced them, i.e. $\mathrm{\Delta }C/\mathrm{\Delta }D={\mathit{\alpha }}_C/{\mathit{\alpha }}_D$. The local endmembers should have been produced before the steady-state achievement by a reaction of similar stoichiometry, but the reaction could have been faster than the observed one or could have occurred in a different place, including in the sediment.

Therefore, in a system defined between only two endmembers, away from atmospheric exchanges, in case of turbulent diffusion mixing, at steady state and with negligible lateral mixing, the reaction-driven approximation allows us to interpret linear variations of TAex versus DICex as a sum of biogeochemical reactions spread all over the water column that can be broken into several reactional stratum. In each stratum, if the local $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ ratio is constant, it corresponds to the apparent stoichiometry of a combination of the biogeochemical reactions occurring in this stratum. In the case of multiple simultaneous reactions in the same stratum, by posing ${\mathit{\alpha }}_C^i$ as the stoichiometry of the ith reaction concerning the reagent C and vⁱ the reaction rate of the ith reaction, we obtain

$${\mathit{\alpha }}_{C}v=\sum _i{\mathit{\alpha }}_C^i{v}^i.$$

To maintain a global reaction rate independent of the species, we have

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }v=\sum _i{v}^i,\\ \text{(4)}& {\displaystyle }\text{so }& {\displaystyle }{\mathit{\alpha }}_C={\displaystyle \frac{\mathrm{1}}{v}}\sum _i{\mathit{\alpha }}_C^i{v}^i.\end{array}$$

Equation (4) indicates that the apparent stoichiometry in a given stratum corresponds to the sum of the stoichiometric coefficients of each reaction weighted by the relative rate of each reaction. Therefore, to estimate the relative rate of each reaction to the observed local changes of TAex, DICex, and AOU (apparent oxygen uptake) or H₂S, a linear combination of reactions is calculated. This combination has to fit three equations (one for each parameter), which allows a maximum of three reactions to be used to solve the system. A limited number of reactions is selected as candidates based on the discussion (see Table 1 and Sects. 4.2 and 4.3). Then, the system is solved with the minimum possible reactions, and the weighted coefficients, $v^i/v$, are calculated.

Table 1 Summary of the main net reactions occurring among the different zones of a redox gradient with reactants starting in equilibrium with the atmosphere, adapted from Soetaert et al. (2007). The calculations assume the Redfield ratio of the organic matter, i.e. γN=0.156 and γP=0.0094. Redox zones describe the different conditions required for completing the reaction.

2.3.2 TA changes indicated by the charge transfer during a reaction

The simplest way to calculate the TA changes induced by an individual reaction is to look at charge transfer induced by the stoichiometry of the given reaction. Indeed, the total alkalinity (TA) corresponds to the quantity of acid added to titrate a solution down to pH 4.5 (Dickson, 1981). It can be described by Eq. (5), with the example of HCl as an acid and B⁻ any titrated base.

$$\begin{array}{}\text{(5)}& \mathrm{HCl}+{\mathrm{B}}^{-}\to \mathrm{BH}+{\mathrm{Cl}}^{-}\end{array}$$

Assuming a complete reaction, the quantity of acid added is equal to the negative charges initially present in the sample consumed plus the positive charges added to the species initially present in the sample. Thus, the total alkalinity corresponds to the loss of negative charges (or gain of positive charges) for species initially present in the sample produced by the pH change from the initial pH = pHini to pH=4.5. Writing $z_i^{\text{pH}}$, the charges held by species i, initially present in the sample at a given pH, we get Eq. (6):

$$\begin{array}{}\text{(6)}& \text{TA}=\sum _i{z}_i^{\text{pH}=\mathrm{4.5}}-\sum _i{z}_i^{\text{pHini}}.\end{array}$$

However, the electroneutrality of water induces

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\sum _i{z}_i^{\text{pHini}}=\mathrm{0}\\ \text{(7)}& {\displaystyle }\text{so }& {\displaystyle }\text{TA}=\sum _i{z}_i^{\text{pH}=\mathrm{4.5}}.\end{array}$$

Equation (7) demonstrates that the total alkalinity is simply the sum of charges that each species present in the sample would have at pH 4.5. From Eq. (7), one can easily deduce the changes of alkalinity from any reaction stoichiometry as soon as the charges of all species at pH=4.5 are known. For example, for a 0.0020 M NaOH dissolution, its TA is 0.0020 M as at pH=4.5, $\text{TA}=\left(\right[{\mathrm{Na}}^{+}]+[{\mathrm{H}}^{+}]-[{\mathrm{OH}}^{-}\left]\right)=\mathrm{0.002}+{\mathrm{10}}^{-\mathrm{4.5}}-{\mathrm{10}}^{-\mathrm{14}+\mathrm{4.5}}=\mathrm{0.0020}$ M. For natural waters, most of the time, the only acid/base species charged at pH=4.5 are ${\mathrm{H}}_{\mathrm{2}}{{\mathrm{PO}}_{\mathrm{4}}}^{-}$ and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$. In that respect, whatever the initial pH and the acid–base equilibrium of species in the sample is, the sum of phosphate species will count negatively, and the sum of ammonium species will count positively. Strictly speaking, at pH=4.5, acid species with pKa between 2.5 and 6.5, such as F⁻ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, would be only partially titrated, and the charge equals their concentration multiplied by a correction factor of $(\mathrm{1}+{\mathrm{10}}^{\text{pKa}-\mathrm{4.5}}{)}^{-\mathrm{1}}$, but this correction can be neglected to a first approximation. Equation (7) is equivalent to those published in Soetaert et al. (2007) or Wolf-Gladrow et al. (2007), whose Eq. (32) can be refined considering that

$$\begin{array}{ll}{\displaystyle }\sum _i& {\displaystyle }{z}_i^{\text{pH}=\mathrm{4.5}}=\left[{\mathrm{Na}}^{+}\right]+\mathrm{2}\left[{\mathrm{Mg}}^{\mathrm{2}+}\right]+\mathrm{2}\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]+\left[{\mathrm{K}}^{+}\right]\\ {\displaystyle }& {\displaystyle }+\mathrm{2}\left[{\mathrm{Sr}}^{\mathrm{2}+}\right]+\mathrm{\cdots }-\left[{\mathrm{Cl}}^{-}\right]-\left[{\mathrm{Br}}^{-}\right]-\left[{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]-\mathrm{\cdots }\\ {\displaystyle }& {\displaystyle }\mathrm{TPO}\mathrm{4}+\mathrm{TNH}\mathrm{3}-\mathrm{2}\mathrm{TSO}\mathrm{4}-\mathrm{THF}-\mathrm{THNO}\mathrm{2}-\mathrm{\cdots }.\end{array}$$

However, Eq. (7) is more concise and more general. For example, in suboxic water, specific species such as polysulfides (as ${{\mathrm{HS}}_{\mathrm{8}}}^{\mathrm{2}-}$, Rickard and Luther, 2007), and in highly productive environments, carboxylic groups from DOC, can be easily added as soon as the bearing charges at pH=4.5 are known.

3 Results

3.1 Water column stratification

High-resolution profiles of salinity and temperature plotted in Fig. 1 show the stratified water column at station 858. Carbonate and redox chemistry are plotted against salinity in Fig. 2. These data result from the 11 CTD casts performed over 1 week during each campaign and correspond to depths ranging from 0.7 to 25 m. Plots against depth generate noisier profiles and are shown in Fig. A2, while plots against salinity follow the water masses. Despite the overall lower salinity due to a near-7-fold river flow increase in 2018 than 2017 (5800 m s⁻¹ in 2018 versus 850 m³ s⁻¹ in 2017), similar zonation of the water column occurred.

Figure 1 Superimposed salinity (in white) and temperature (in black) profiles over the 11 casts for each campaign.

Figure 2 Superimposed carbonate and redox chemistry profiles of 11 casts done in August 2017 and in August 2018 at station 858 of the Chesapeake Bay.

A surface layer sampled only in 2018 is visible in the top 3 m depth (Fig. 1). It presents highly variable temperature and oxygen concentration and oversaturation of pCO₂ (Fig. 2), indicating emission to the atmosphere. Below, at 3 m depth, a subsurface layer (named primary production zone or PP in Fig. 2) is characterised by a high amount of O₂ (about or above 100 % saturation), high pH (about 8; 8.11 ± 0.07, n=13 in 2017 and 7.94 ± 0.08, n=14 in 2018) and high day-to-day temperature variation (above 1 ^∘C between different days). The layer presents relatively low pCO₂ (505 ±75 µatm, n=13 in 2017 and 770 ± 130 µatm, n=14 in 2018), with minimal values at 110 µatm in 2017 and 205 µatm in 2018 which are below the atmospheric pCO₂ of 407 µatm (Chen et al., 2020). This signature corresponds to primary production (PP). Fluorescence (not shown) correlates with pH as expected for primary production ($\text{pH}=\text{Fluo}\left(\mathrm{mV}\right)\times \mathrm{13}+\mathrm{7.14}$, r²=0.8 in 2017 and r²=0.9 in 2018). DICex and TAex reach their minimal value in the surface layer. While DICex minimum is similar between the two campaigns (78 ± 17 µM, n=13 in 2017 and 76 ± 10 µM, n=14 in 2018); TAex minimum is much lower in 2017 (7 ± 3 µM, n=12) than in 2018 (48 ± 7 µM, n=14).

Below, with increasing salinity, an important increase of pCO₂ accompanying the decrease of O₂, pH and temperature is visible. A relatively invariable low oxygen zone (called ILO in Fig. 2) is here defined by the salinity invariance of O₂ concentrations at about 24 ± 3 µM (n=15) in 2017 and 105 ± 5 µM (n=16) in 2018. Other species are also relatively stable for this depth, such as pCO₂ at about 3250 ± 150 µatm (n=10) in 2017 and 2590 ± 100 µatm (n=16) in 2018 and pH at about 7.28 ± 0.01 (n=13) in 2017 and 7.37 ± 0.02 (n=16) in 2018. Deeper, where the oxygen is not detectable (≲1 µM) and in absence of free sulfide, the so-called anoxic zone corresponds to a pH minimum at 7.24 ± 0.01 (n=40), similar in 2017 and 2018, that generates a pCO₂ maximum. The deepest layer is a sulfidic layer, [H₂S]=11.2 ± 2.8 µM (n=36), in which the pH seems quite stable at 7.33 ± 0.01 (n=29) and 7.32 ± 0.01 (n=9) in 2017 and 2018 respectively. The Ca²⁺ concentrations observed by Su et al. (2021) and during the 2018 cruise (data not shown) vary linearly with salinity (calcium excess stays below 200 µM or 10 % of total Ca). Assuming similar behaviour in 2017, calculations show that the whole water column (except four samples from the PP zone) is under saturated ($\mathrm{0.36}<\mathrm{\Omega }\text{cal}<\mathrm{1}$; mean=0.68) with respect to calcite in 2018, while undersaturation is only valid below S=10 in 2017. The main changes between the two campaigns correspond to a greater oxygen penetration in 2018, especially visible in the ILO zone. Additionally, a surface layer (with salinity below 3) is visible above the primary production zone in 2018 but is related to the more superficial sampling in 2018.

3.2 Intermediate redox species

The development of the anoxic zone during summer (Su et al., 2021) and the regularity of this development over the years (Sholkovitz et al., 1992; Trouwborst et al., 2006; Lewis et al., 2007; Cai et al., 2017; Oldham et al., 2017a) requires the presence of species able to rapidly oxidise the H₂S mixing upward and to reduce the O₂ mixing downward. The three main redox couples known to play this role, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}/{{\mathrm{NO}}_{\mathrm{2}}}^{-}$, ${\mathrm{MnO}}_x/{\mathrm{Mn}}^{\mathrm{2}+}$ and ${\mathrm{Fe}}^{\mathrm{3}+}/{\mathrm{Fe}}^{\mathrm{2}+}$, are described in Fig. 3 by the superimposition of all cast results against salinity. Four representative casts are plotted in Fig. A3. The primary production zone is depleted in dissolved Mn(II) and Mn(III) and shows an average value for solid MnOx (0.7 ± 0.1 µM, n=4) and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (2.1 ± 0.1 µM, n=18). In 2018, iron speciation was investigated in the primary production zone, showing an important iron pool (1.6 ± 0.1 µM, n=9) dominated by solid Fe(II) (95 % ± 2 %). Below, in the ILO zone, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reaches a maximum plateau at 5.3 ± 0.2 µM, n=11 in 2017 that is not visible in 2018; MnOx dominates the Mn pool in the ILO zone with a concentration of 1.8 ± 0.1 µM, n=14 in 2017 and 1.3 ± 0.2 µM, and n=3 in 2018. Fe(II) represents only 49 %± 13 %, n=4 of the iron pool (for a total iron concentration of 2.0 ± 0.4 µM, n=4) in 2017, while it represents 87 % ± 5 %, n=4 of the iron pool in 2018 (for a total concentration of 1.2 ± 0.2 µM, n=4). Just below oxygen depletion, in the anoxic zone, MnOx reaches a maximum (2.0 ± 0.1 µM, n=25 in 2017 and 2.9 ± 0.4 µM, n=12 in 2018). During both campaigns, the sulfidic layer is characterised by the absence of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and MnOx, while dissolved manganese concentration increases up to 7.7 ± 0.8 µM, n=15, and the iron pool increases to 5.0 ± 0.7 µM, n=13. No Mn(III) was detected with the porphyrin kinetics method (Thibault de Chanvalon and Luther, 2019), but about 30 % of the total dissolved manganese flocculated after acidification down to pH 1.5 when analysed in 2017, indicating the existence of Mn(III).

Figure 3 Mn, Fe and nitrite profiles built by superimposition of 11 casts for each campaign.

4 Discussion

In terms of salinity, station 858 of the Chesapeake Bay shows very similar water column zonation between summer 2017 and summer 2018, despite a 10-fold difference in freshwater discharge rates (Figs. 1 and 2). Major features are first, a surface layer characterised by intense atmospheric exchange that was only sampled in 2018. Below, at about 3 m depth, there is a subsurface layer associated with high primary production (PP zone) with high pH of about 8, oversaturation of dioxygen and low CO₂ partial pressure (down to 110 µatm). Below, a low oxygen layer with invariant concentration of most species surveyed (the invariant low oxygen (ILO) zone) is characterised by significant nitrite accumulation in 2017 (Fig. 2), probably due to oxidation of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ diffusing upward and/or produced by in situ remineralisation. This feature is not visible in 2018 probably because the higher O₂ concentration in 2018 accelerates nitrite oxidation into nitrate and prevents any significant accumulation. Below, the anoxic zone, with neither O₂ nor H₂S detectable, is characterised by an increase of MnOx concentration and a pH minimum. This MnOx maximum can be explained by the upward diffusing Mn²⁺ that is biologically oxidised by the downward diffusing O₂, though at low, undetectable, concentration (Clement et al., 2009). Additionally, Mn²⁺ could be oxidised by the nitrite or the nitrate (not measured) diffusing downward (thermodynamically favourable, Luther, 2010). Compared to the ILO zone, the anoxic MnOx maximum corresponds to an increase of 0.2 µM in 2017, while it is much more marked in 2018 with an increase of 1.6 µM (Fig. 3). This difference could come from a faster Mn²⁺ oxidation produced by the steeper oxygen gradient above the anoxic zone and by the thinness of the anoxic layer in 2018. Finally, in the deeper sulfidic layer, the MnOx disappearance corresponds to the Mn²⁺ increases (Fig. 3) according to the reduction of settling MnOx by H₂S. The concentration increase of the manganese pool and of the iron pool with depth in the anoxic and sulfidic layers probably results from important sedimentary efflux.

4.1 Validity of the reaction-driven approximation

The simultaneous and high-resolution sampling of multiple carbonate parameters and redox species gives us the rare opportunity to investigate in detail the interaction between carbonate species and redox-sensitive elements. In the Chesapeake Bay, main changes of DIC and TA can be explained by mixing between upstream and oceanic endmembers (Appendix, Fig. A1). This “endmember-driven” interpretation leads to the calculation of an excess of DIC and TA, DICex and TAex, relative to the mixing line, as shown in Fig. 4a. At station 858, the steep gradient observed, for example the pH and pCO₂ gradients in the PP zone; the O₂ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ gradients above the anoxic zone; and the Mn, Fe and H₂S gradients at depth, suggest that ongoing in situ processes control the changes of concentrations and dominate the time-dependent endmember variability or the mixing with an unknown third endmember. Additionally, the TAex versus DICex plot (Fig. 4a) shows steep changes of direction, such as in the case of straight lines between local endmembers maintained by ongoing reactions, while a preponderance of mixing would produce more progressive changes. Finally, at depth, Fig. 4a shows a similar slope for both years studied ($\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}=\mathrm{2.4}$) rather than a similar TAex and DICex concentration, which reinforces the validity of the reaction-driven approximation.

Figure 4 Description of the TAex$/$DICex$/$AOU or H₂S system. Panels (a), (b) and (d) show samples measured along with interpretative slopes discussed in the main text. Panel (c) shows the theoretical slope of TAex$/$DICex from combination of reactions presented in Table 1.

Assuming (1) that mixing is efficiently described by vertical turbulent diffusion mixing; (2) that the measured concentrations correspond to a steady state – no changes observed over the 1 week sampling; (3) that no additional endmember contributes significantly to the excess calculation, in particular that the samples are isolated from atmospheric exchanges; and (4) that lateral mixing is negligible, which is equivalent to the lateral invariance of the system – as in the stratified water column of station 858, the reaction-driven approximation (Sect. 2.3.1) permits interpretation of the concentration changes as a linear combination of the stoichiometry of several chemical reactions (Eq. 5). The weighted coefficients of each reaction are equal to the rate of each reaction relative to the sum of the rate of all occurring reactions, $v^i/v$ (Eq. 4). Accordingly, this interpretation does not identify reactions with minor impact on the carbonate cycle (Eq. 7) or reactions that cancelled each other later during the journey. For example, ΔTAex introduced by PP is frequently cancelled by a similar amount or even higher ΔTAex by AR. The starting point corresponds to the upstream (S=1.5) endmember at DICex=0 and TAex=0. Figure 4a shows a slight DICex enrichment at TAex=0, which reflects a bias from the endmember calculation.

4.2 Identification of major reactions in oxygenated water column using $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ signature

In oxygenated water, the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ analysis can be supported and strengthened by looking also at the $\mathrm{\Delta }\text{AOU}/\mathrm{\Delta }\text{DICex}$ ratio (Fig. 4b) and establishing a $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ signature for each water layer. Then, a linear combination of up to three reactions can be fit to the observed $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ signature. The main candidate reactions to be combined are aerobic respiration (AR), primary production (whose overall mass balance equation is here summarised as negative AR), carbonate dissolution (CD) and nitrification of ammonium (Nit) as presented in Table 1 and in Fig. 4c. In 2017, the whole oxygenated zone, including the PP zone and the ILO zone, is characterised by a $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=\mathrm{0.1}/\mathrm{1}/\mathrm{1}$ (Fig. 4a and b), which corresponds to the occurrence of only net aerobic respiration (AR) (theoretical values are $\mathrm{0.15}/\mathrm{1}/\mathrm{1}$, see Table 1 and Fig. 4c). Note (a) that “net aerobic respiration” indicates that primary production is possible at a significant rate but slower than AR and (b) that this $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ signature indicates weak nitrification following respiration. In case of full nitrification, the theoretical slopes should be $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=-\mathrm{0.167}/\mathrm{1}/\mathrm{1.31}$, as proposed by Zeebe and Wolf-Gladrow (2001). A combination of 1 AR and 0.025 Nit (nitrification of 16 % of the produced NH₃) improves the modelled value to $\mathrm{0.1}/\mathrm{1}/\mathrm{1.05}$. The relatively slow nitrification can be explained by slow kinetics of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ oxidation, with a half-life time estimated between a few days in estuaries (Horrigan et al., 1990) to multiple years in coastal environments (Heiss and Fulweiler, 2016). As a comparison, other ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ fates, such as adsorption, lead to an ammonium half-life time of about a few minutes (Alshameri et al., 2018) to a few hours (Raaphorst and Malschaert, 1996) depending on the concentration of fine particles. Additionally, algae are known to use ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ as an N source (Raven et al., 1992), and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ can be directly assimilated by heterotrophic organisms.

In 2018, fresh water masses brought by the exceptional flood modified the carbonate system equilibrium. First, a low-salinity layer with pCO₂ at 1540 µatm overlays the primary production layer (Fig. 2), preventing the uptake of atmospheric CO₂ by primary production as was observed in 2016 (Chen et al., 2020). Just below the surface layer, in the PP zone, the lock down of atmospheric exchange by the low-salinity layer produces supersaturation of trapped O₂ (Fig. 2). In Fig. 4a and b, this process translates into a vertical distribution at DICex=40 µM with $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=\mathrm{1.37}/\mathrm{0}/-\mathrm{1}$. This original signature can be modelled by the combination of simultaneous carbonate dissolution (CD), the water column being undersaturated, and PP, no important turbidity was visible as modelled by Cerco et al. (2013), in equal proportion (second line in Table 2); the carbonate dissolution buffers the DIC consumption by the PP. Note that the ratio between $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{AOU}$ implicates an important nitrate assimilation superior or equal to the amount of N required for the PP, as modelised by negative nitrification in Table 2.

Table 2 Linear combination of reactions from Table 1 that fit the observations (see text for details, H₂O molecules are omitted).

In 2018, below the PP zone down to DICex=240 µM, the beginning of the ILO zone, the TAex increases significantly with a $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=\mathrm{0.67}/\mathrm{1}/\mathrm{0.67}$ (Fig. 4a and b), incompatible with AR. To explain this signature, a contribution of CD superimposed on AR seems most likely. A linear combination fitting leads to 0.5 CD for 1 AR (third line in Table 2) and results in $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=\mathrm{0.77}/\mathrm{1}/\mathrm{0.67}$. An explanation of its occurrence solely in 2018 could be the increase in carbonate-rich suspended material at high flow conditions (Su et al., 2021). Excess of Ca²⁺ compared to the mixing line with oceanic end member (not shown) indicates that up to 200 µM of Ca²⁺ is produced in the oxygenated layer. Deeper, in the ILO zone, $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=\mathrm{0.2}/\mathrm{1}/\mathrm{1.25}$ (Fig. 4a and b), which results mainly from AR ($\mathrm{0.15}/\mathrm{1}/\mathrm{1}$) with the possible addition of CD and Nit, the exact signature being fitted for 0.54 CD and 0.46 Nit for 1 AR (fourth line in Table 2), as a continuum of AR and CD relative rates from the overlaying layer. This important nitrification is also in good agreement with the nitrification observed in the PP zone during this campaign, the lack of nitrite build-up in the ILO zone and the relatively high oxygen concentration (at 105 µM) in the ILO zone in 2018 which is able to sustain nitrification. Overall, our results show that the higher river flow of 2018 increases carbonate dissolution for the top 5–10 m water depth that is superimposed on primary production or aerobic respiration. Additionally, in 2017 the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ system indicates weak nitrification, while in 2018 significant nitrification in the ILO and PP zones is suggested by the reaction-driven approximation. The role of nitrification in explaining TAex depletion is only hypothetical since no direct measurements of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ were performed. In particular, TAex depletion is particularly intense during high flow, high suspended particles season and could be produced by ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ adsorption to the particles rather than by nitrification.

4.3 Identification of major reactions in the anoxic water column using the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ signature

In the absence of oxygen, the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ will be used for rates calculation. In the anoxic zone, the signature is similar in summer 2017 and 2018 at $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}=\mathrm{2.4}/\mathrm{1}/\mathrm{0}$. The reaction-driven approximation of this signature is more difficult than in oxidised water because more reactions are known to occur simultaneously in the absence of oxygen. However, to fit with the reaction-driven approximation, it is not necessary to describe each reaction step but only the overall changes concerning the journey of a water mass over different redox conditions resulting in a combined result or reaction. This approach has recently been proposed for FeS burial by Hiscock and Millero (2006), Rassmann et al. (2020) or Su et al. (2020b). A scenario combining sulfate reduction (SR) is particularly attractive since SR represents the main carbon remineralisation pathway in absence of oxygen. However, a combination of SR with CD would result in a $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ between 1.15 and 2 (see Table 1 or Fig. 4c) and fails to reach the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ of 2.4. Moreover, SR alone underestimates the importance of the H₂S oxidation pathway that can consume all the alkalinity produced during SR. For example, SR follow by oxygenated oxidation results in a $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ signature equal to AR only. In the Chesapeake Bay, H₂S oxidation is critical since no H₂S is measurable in the anoxic zone, while the gradient at the sediment–water interface indicates high H₂S sedimentary efflux (Fig. 2).

Generalising these observations, recent efforts to build an alkalinity budget on the global scale (Hu and Cai, 2011; Middelburg et al., 2020) highlight that the alkalinity produced by anaerobic respiration corresponds to the uncharged species produced, mostly in solid or gaseous phases. Indeed, the alkalinity changes produced during a natural reaction equal the “charge transfer” from species having some charge at pH=4.5, such as ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$, to a species that would lose its charges at pH=4.5, mainly ${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$, that is not counted in the alkalinity calculation (see Eq. 12). Although correct, this approach tends to neglect the roles of Fe and Mn oxides (Middelburg et al., 2020) since their transformation from (oxyhydr)oxides into sulfur or carbonate solid species does not involve any charge transfer. When looking in detail at these processes, the metal oxides are critical since they are the main H₂S oxidation pathway that does not regenerate H₂SO₄ but rather produces S⁰, ${{\mathrm{S}}_n}^{\mathrm{2}-}$ or FeS (Findlay et al., 2014; Avetisyan et al., 2021), which limits alkalinity consumption.

To build a pool of candidate reactions for the fitting, first, dissolved species at too low a concentration (e.g. Mn²⁺_aq, Fe²⁺_aq) to be a net reagent to affect the carbon cycle at steady state are not taken into account. These species are usually recycled rapidly and hold the role of catalyser or electron shuttle between other redox species and did not reach 10 µM during the campaigns (Fig. 3). Second, many minerals are expected to be at low concentration or thermodynamically not favoured, and their associated reactions are neglected (e.g. iron phosphate, ferrous or manganous oxide, metal sulfur clusters, MnS, FeCO₃, adsorption processes, reverse weathering). Therefore, only aqueous species with important stock concentrations (that can exceed 0.1 mM in anoxic water) are taken into account, i.e. ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$, Ca²⁺, H₂S and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ together with gaseous (N₂, CO₂) and main solid phases (FeS₂, FeS, S⁰, MnCO₃, FeOOH, MnO₂). Third, the combination of many carbon remineralisation reactions with re-oxidation reactions or the net result produces a net chemical equation equal to another remineralisation reaction. As an example, the chemical equation of SR followed by H₂S oxidation with oxygen is equal to the equation of aerobic respiration: the proposed model confounds both pathways because the resulting chemical changes are similar.

Table 1 lists the resulting combined reactions, and the calculated $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ slopes are represented in Fig. 4c. For the anoxic zone, the nitrate respiration can be associated with N₂ production (NR-N2) or with NH₃ production (NR-NH₃), and the manganese oxide respiration produces carbonate precipitation (MnR-MnC). For the sulfidic zone, SR can occur alone, producing a build-up of H₂S. However, at a certain point, H₂S gets significantly oxidised either by MnO₂ which produces MnCO₃ and S₀ (SR-SMnC) or by FeOOH, producing FeS and S₀ (SR-SFeS) and ultimately FeS₂ (SR-FeS2). Direct respiration of FeOOH is also taken into account, but as the only final Fe product in the model is FeS or FeS₂, it has to be accompanied by some SR (FeSR-FeS).

Figure 4c demonstrates that the slope of $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}=\mathrm{2.4}$ can be obtained for any reaction in combination with MnR-MnC. Combinations without MnR-MnC, however, lead to a negative SR, whose overall equation could be interpreted as a possible small participation of anoxygenic phototrophic (purple) bacteria (Findlay et al., 2015, 2017), but are not considered further, as the amount of ΔTAex involved would be tiny. Therefore, in the absence of nitrate, oxygen and H₂S, only a combination of MnR-MnC with SR-SMnC (producing S⁰, fifth line in Table 2), SR-FeS (producing FeS, sixth line in Table 2) or CD (releasing Ca²⁺, seventh line in Table 2) gives the particularly high $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ of 2.4. S₀ was not measured during our campaign, but it has been previously reported in this water column (Findlay et al., 2014), and the S₀ produced by SR-SMnC can react with FeS to form FeS₂. All the three identified combinations require a critical role of MnO₂. Since Fig. 3 does not indicate any clear reaction for iron while the steep gradient of MnOx in proximity to the H₂S rich layer suggests reaction between MnOx and H₂S, the combination with SR-SMnC is the most likely (fifth line in Table 2). Deeper, the vertical gradient of sulfide suggests that part of the H₂S came by diffusion from the sediment's porewater (Fig. 2). The assumptions required for the reaction-driven approximation are still valid as soon as steady state is maintained by ongoing reactions, even if one of the local endmembers has not been sampled, since it is probably located in the sediment. In the presence of sulfide, the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ signature is $\mathrm{2.4}/\mathrm{1}/\mathrm{1.2}$ in 2017 and $\mathrm{2.4}/\mathrm{1}/\mathrm{3.2}$ in 2018 (Fig. 4a and d) and can be explained by the same combination of reactions without complete oxidation of H₂S from SR to take into account the build-up of H₂S (Table 2).

4.4 Comparison with other studies

The high observed ratio of $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}=\mathrm{2.4}$ seems very specific to the Chesapeake Bay. Moreover, the reaction-driven interpretation can be applied to other published datasets for which the $\mathrm{\Delta }\text{TA}/\mathrm{\Delta }\text{DIC}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ system can be calculated (Table 3). In the water column, most of the available datasets are not suitable for the reaction-driven approximation, since either they focus on surface water where DIC and TA are strongly impacted by atmospheric exchange or the water masses change too fast to consider that reactions dominate over water mixing. However, porewater measured in the Gulf of Mexico has $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}=\mathrm{1.15}/\mathrm{1}/\mathrm{0.53}$ (Hu et al., 2010), as expected when sulfate reduction is associated with H₂S accumulation (SR reaction; $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}=\mathrm{1.15}/\mathrm{1}/\mathrm{0.5}$; Table 1). Similarly, Hiscock and Millero (2006) report $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}=\mathrm{1.3}/\mathrm{1}/\mathrm{0.5}$ in the western Black Sea close to the SR signature. In the Baltic Sea sediment, $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}=\mathrm{1.3}/\mathrm{1}/\mathrm{0.07}$ was reported (Lukawska-Matuszewska, 2016), which is close to the expected signature in the case of important H₂S consumption by Fe oxides and consequent precipitation as pyrite (SR-FeS2 reaction; $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}=\mathrm{1.2}/\mathrm{1}/\mathrm{0}$). In Rhône river prodelta sediments, the reported $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ is $\mathrm{1}/\mathrm{1}/\mathrm{0}$, which can be related to the $\mathrm{1}/\mathrm{1}/\mathrm{0}$ signature of FeSR-FeS (Table 1) that is expected in iron-rich sediments, with high sedimentation rates preventing FeS₂ formation in the pore water (Rassmann et al., 2020). In permeable, carbonate-rich sediments, the reported signature of $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ in Hawaii sands is $\mathrm{0.86}/\mathrm{1}/\text{ND}$ (Drupp et al., 2016). The lack of salinity and oxygen datasets prevents further model fits, but the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ is below 2.4. For oxygen-depleted data from the whole Chesapeake Bay described in Su et al. (2020a), the signatures are $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}=\mathrm{0.2}/\mathrm{1}/\mathrm{1}$ in the presence of oxygen (typical signature of AR) and $\mathrm{0.8}/\mathrm{1}/\text{ND}$ in the absence of oxygen that could correspond to NR-N2 or a combination of CD and Nit. Overall, this bibliographic survey highlights the effectiveness of the reaction-driven approximation to identify preponderant reactions controlling the carbon cycle and puts into perspective the originality of the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}/\mathrm{\Delta }\text{AOU}$ signature of $\mathrm{2.4}/\mathrm{1}/\mathrm{0}$ observed in the Chesapeake Bay.

Table 3 Overview of the $\mathrm{\Delta }\text{TA}/\mathrm{\Delta }\text{DIC}/\mathrm{\Delta }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ signature observed in different environments.

^∗ Lack of data on oxygen concentration prevent any interpretation. ${}^{\ast \ast }$ Important air–water exchange prevents any reaction-driven approximation.

4.5 Local budget

While the reaction-driven approximation indicates a dominant role of the SR-SMnC reaction, this possibility needs to be validated by looking at the saturation state of rhodochrosite (the main MnCO₃ mineral) and looking at the mass budget between MnOx consumed and TAex produced. The rhodochrosite saturation (Luo and Millero, 2003) is always below 0.3 in our samples, which stands against the occurrence of in situ SR-SMnC reaction. When inspecting the mass budget, the 88 to 155 µM of MnO₂ required to produce the 100 µM TAex increase observed (Fig. 4a) is 1 order of magnitude higher than the observed MnOx or Mn²⁺ concentration (Fig. 3). This mass budget discrepancy cannot be solved by invoking suspended material, since the 88 µM of MnO₂ would require a suspended material concentration of about 4.4 g L⁻¹ (assuming an average concentration of 20 µmol g⁻¹ of Mn), which is again 1 or 2 orders of magnitude higher than the 0.01–0.1 g L⁻¹ usually found in the Chesapeake Bay (Cerco et al., 2013). However, a fast settling rate could satisfy and explain the discrepancy between water and solid concentration. But another process dephasing aqueous from solid reaction products is also possible at station 858, since the dissolved phase could have moved up, rather than the particles settling down. In this case, the SR-SMnC reaction was not happening only in the water column of the Chesapeake Bay, and part of the TAex and DICex pool could have been produced in the sediment during the previous year (e.g. Aller, 2014), then diffused out of the sediment simultaneously with other reduced elements as the summer begins. Indeed, many previous studies at station 858 (e.g. Sholkovitz et al., 1992) explained the seasonality of anoxia with an upward move of the redox front from the sediments to bottom waters during the start of summer. Important sedimentary efflux of H₂S, Fe_aq and Mn_aq were still visible during both of our August campaigns. Therefore, the 100 µM TA pool does not fit with the ambient Mn²⁺ or MnO₂ in the water column but rather with the MnCO₃ deposited in the sediment. The sedimentary solid Mn stock of the Chesapeake Bay is particularly important, up to 70 µmol g⁻¹ at station 858 (Sinex and Helz, 1981) compared to an average value of 15 µmol g⁻¹ for the upper continental crust (Rudnick and Gao, 2003). Indeed semi-enclosed basins are known to concentrate manganese at the deeper sediment (Thamdrup and Dalsgaard, 2000; Lenstra et al., 2020). Recent investigations at a close station (ET 5.1; $\mathrm{38}{}^{\circ }{\mathrm{48.36}}^{\prime }$ N; $\mathrm{75}{}^{\circ }{\mathrm{54.66}}^{\prime }$ W) in the Chesapeake Bay (Lenstra et al., 2021) report that about 60 % of the surface sedimentary Mn pool is MnOx (acid ascorbic extractable) and 25 % is Mn carbonate (1 M HCl extraction). Assuming a porosity of 0.8 and a bulk solid density of 2.6, the sedimentary pool corresponds to 35 mM of manganese, which largely exceeds the 88 µM required to produce the 100 µM TAex increase. Therefore, the Chesapeake Bay sediment is particularly rich in manganese and could host important SR-SMnC reactions in the surface pore water, whose soluble products diffuse up to the water column during summer and could bear with them the high $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ signature observed.

The assumptions detailed in Sect. 4.1 permit the reaction-driven approach to be reconsidered in regard to this sediment efflux scenario. Indeed, the sediment efflux does not need to be considered as an additional endmember; that would violate the third assumption, since its salinity and the pore water concentrations results from the upstream and oceanic endmembers superimposed to geochemical reactions. In Sect. 2.3.1, we point out that Eq. (1) was valid in between each local endmember and that the straight lines on the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ plot between them indicate that the local endmembers are maintained in steady state by a chemical reaction with a similar stoichiometry. These results indicate that, even if most of the MnCO₃ was produced when the local endmembers were localised in the sediment, (a) their migration does not alter the TAex$/$DICex signature, and (b) the chemical reaction that produced them is still ongoing at a sufficient rate to maintain a steady state characterised by the steep changes of direction observed in Fig. 4a.

5 Conclusion

The reaction-driven approximation proposed in this study is a powerful interpretative framework that can identify the major reactions controlling the carbonate cycle. In the Chesapeake Bay, similar redox stratification can support varying intensity of carbonate dissolution, absent as in 2017 or important as in 2018. In 2018, the reaction-driven approximation also suggests an important role of the N cycle to consume alkalinity, while in 2017, nitrification was limited to 16 % of the nitrogen mineralised. The summer anoxia observed in the Chesapeake Bay is characterised by an exceptionally high $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ of 2.4 which has never been reported in anoxic water columns or sediment pore waters. The reaction-driven approximation suggests that it comes from sulfate reduction with almost complete hydrogen sulfide oxidation by MnO_x followed by MnCO₃ precipitation. This interpretation is supported by the important manganese dynamics observed in this and previous papers (Oldham et al., 2017a; Sholkovitz et al., 1992; Trouwborst et al., 2006). However, most of the reactions would have occurred in the sediment from which the components of the $\mathrm{\Delta }\text{TAex}/\mathrm{\Delta }\text{DICex}$ of 2.4 diffused into the bottom water with the redox front during the set-up of summer water column anoxia. Our study demonstrates that the manganese cycle can have a strong impact on alkalinity, as it prevents H₂S oxidation to ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ and favours sulfur burial.

Appendix A

A1 Demonstration of Eq. (1)

First, Eq. (A1) describes the behaviour of a solute, in the case of turbulent diffusion mixing, with D_s the effective diffusion coefficient, superimposed on a chemical reaction of rate v as described in Eq. (A2).

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=D_s{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}C}{\mathrm{d}{x}^{\mathrm{2}}}}+{\mathit{\alpha }}_{C}v\text{(A2)}& {\displaystyle }{\mathit{\alpha }}_{C}C+{\mathit{\alpha }}_{D}D={\mathit{\alpha }}_{E}E+{\mathit{\alpha }}_{F}F\end{array}$$

Assuming steady state, i.e. $\mathrm{d}C/\mathrm{d}t=\mathrm{0}$, and applying Eq. (1) to two different species sharing the same reactions at rate v on each point and α_D≠0, we can express the rate of the reaction as

$$\begin{array}{}\text{(A3)}& -{\displaystyle \frac{v{\mathit{\alpha }}_C}{D_s}}={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}C}{\mathrm{d}{x}^{\mathrm{2}}}}={\displaystyle \frac{{\mathit{\alpha }}_C}{{\mathit{\alpha }}_D}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}D}{\mathrm{d}{x}^{\mathrm{2}}}}.\end{array}$$

Equation (A3) is also true for any linear combination of a solute with a conservative element such as the salinity, S, because α_S=0, hence (with k any real)

$${\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}C+kS}{\mathrm{d}{x}^{\mathrm{2}}}}={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}C}{\mathrm{d}{x}^{\mathrm{2}}}}+k{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}S}{\mathrm{d}{x}^{\mathrm{2}}}}={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}C}{\mathrm{d}{x}^{\mathrm{2}}}}.$$

By integration of Eq. (A3) with respect to x on a portion of space where ${\mathit{\alpha }}_C/{\mathit{\alpha }}_D$ is constant, it becomes Eq. (1).

A2 Validity of Eq. (2) when the two endmembers have similar concentrations

From Eq. (A3), and assuming that the rate v can be described by a polynomial function of x, with a_i, any real coefficient, we obtain

$${\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}C}{\mathrm{d}{x}^{\mathrm{2}}}}=-{\displaystyle \frac{{\mathit{\alpha }}_C}{D_S}}\sum _i{a}_i{x}^i.$$

That can be solved analytically and gives the following equation:

$${\displaystyle \frac{\mathrm{d}C}{\mathrm{d}x}}={\displaystyle \frac{{\mathit{\alpha }}_C}{D_S}}\sum _i\left[{\displaystyle \frac{a_i}{(i+\mathrm{1})(i+\mathrm{2})}}(\mathrm{1}-(i+\mathrm{2}\left)x^{i+\mathrm{1}}\right)\right]+\mathrm{\Delta }C.$$

Here, ΔC is the difference of concentration between the two endmembers. With similar reasoning on D and with $\mathrm{\Delta }C=\mathrm{\Delta }D=\mathrm{0}$, we find Eq. (1) with G=0.

A3 Additional figures

Figure A1 Variation of total alkalinity (TA) during oceanic and river mixing. Dashed lines represent the theoretical DIC and TA if only mixing occurs.

Figure A2 Superimposed carbonate and redox chemistry profiles over 11 casts against depth.

Figure A3 Examples of redox and metal chemistry profiles. Note that the Fe(III)_s and Fe(III)_aq scales are zoomed in by 3 and the O₂ scale is zoomed out by 100.