Sulfate enrichment in estuaries of the northwestern Gulf of Mexico: The potential effect of sulfide oxidation on carbonate chemistry under a changing climate

Water quality parameters from 2000 to 2020 were used to identify the spatial and temporal sulfate variations in estuaries of the northwestern Gulf of Mexico. Sulfate enrichment relative to conservative mixing was found to be associated with a low river discharge period from 2012 to 2014 in all estuaries. Based on reaction stoichiometry, sedimentary sulfide oxidation holds significant potential for reducing the alkalinity in estuarine waters. However, during this extreme drought, alkalinity enrichment was also occasionally observed in some of the southern estuaries along with sulfate enrichment, and when alkalinity depletion occurred, the magnitude of depletion was usually much less than what would be expected based on sulfide oxidation alone. This discrepancy can be partially explained by carbonate dissolution and other proton removal pathways (e.g., Fe‐oxide dissolution), and by the uncertainties in the model used to estimate alkalinity enrichment/depletion. Under a changing climate, the close coupling between river discharge variation and estuarine sulfate dynamics will significantly impact estuarine carbonate chemistry.

As atmospheric CO 2 increases due to anthropogenic activities, oceanic uptake of atmospheric CO 2 has lowered seawater pH, altered ocean biogeochemistry, threatened the health of aquatic organisms, and impacted fisheries resources and economies.On a local scale, ocean acidification also affects the carbonate chemistry of estuaries by transporting low pH waters into estuaries (Feely et al. 2010;Hu and Cai 2013).Due to the strong horizontal stratification, long water residence time, eutrophication, and a weak acid-base buffer capacity, the CO 2 -induced acidification effect is magnified in many estuaries (e.g., Scheldt Estuary and Chesapeake Bay) (Hofmann et al. 2009;Cai et al. 2017).
In addition to river-ocean mixing, biogeochemical processes including but not limited to primary production/ remineralization, nitrification/denitrification, iron reduction/ oxidation, and carbonate dissolution/precipitation all contribute to alkalinity variations in an estuary (Hunt et al. 2022).Sulfate reduction and sulfide oxidation are also important reactions that contribute to coastal and estuarine alkalinity balance (Hu and Cai 2011), even though sulfate increase can also be caused by river discharge with high sulfate concentrations (Kwiecinski 1965).A temporary excess of sulfide oxidation over sulfate reduction enriches sulfate, and sulfate enrichments have been observed in tidal creeks and upper layers of estuarine and marsh sediments (Luther et al. 1982;Matson 1982;Lord and Church 1983).In addition, estuarine water sulfate enrichment and removal fluctuations have been reported in Pamlico and Neuse River estuaries, North Carolina, USA (Matson and Brinson 1985), and a potential impact on water column carbonate chemistry could be inferred, although no carbonate chemistry parameter was determined.Recently, a lab-based incubation study also demonstrated that sulfate production along with alkalinity consumption occurred in the shallow, semiarid Mission-Aransas Estuary located in the northwestern Gulf of Mexico during droughts (Dias et al. 2023).
Sulfate enrichment and depletion can be calculated by deviations from the conservative sulfate-salinity relationship.These deviations do not necessarily persist due to rapid sulfur cycling or hydrological and tidal exchange (e.g., Luther et al. 1982).Estuaries in the northwestern Gulf of Mexico are characterized by latitudinal gradients in both rainfall and freshwater inflow, which result in a southward increasing water residence time (from less than a month to over a year), increasing average salinity, and decreasing salinity gradient within the estuaries, even though all the estuaries share similar physiogeographical characteristics (Montagna et al. 2012).These distinct characteristics may lead to different biogeochemical responses to climate variations.
To understand sulfur cycling and its potential effect on estuarine carbonate chemistry under the changing climate, we examined the spatial and temporal relationships between river discharge, estuarine water sulfate, and alkalinity in the northwestern Gulf of Mexico estuaries using data collected from 2000 to 2020.We interpreted the variabilities of these solutes based on biogeochemical processes and reaction stoichiometry.We found that sulfate enrichment (i.e., positive deviation from river-ocean conservative mixing) generally coincided with a low river discharge (drought) period (2012)(2013)(2014), indicating significant impacts on estuarine carbonate chemistry.Alkalinity depletion occurred from 2012 to 2014 in most cases, although its loss was much smaller than what would be expected based on the stoichiometry of the sulfide oxidation reaction alone.Part of the discrepancy could be explained by sediment carbonate dissolution, which acted as a sink for the acid produced.Besides the large uncertainties in the river endmember alkalinity, freshwater evaporation also regulated the nonconservative behavior of estuarine water alkalinity, especially during an extreme drought in the southern estuaries.Climate models project an unprecedented drought risk in the southwestern USA during the latter half of the 21 st century (Cook et al. 2015).With drier climate conditions in the future, this study suggests an increasing sulfate enrichment trend caused by the oxidation of reduced sulfur buried in sediments, which may be accompanied by the dissolution of recently deposited sedimentary carbonate along with water alkalinity loss in semiarid estuaries.

Study area
Most northwestern Gulf of Mexico estuaries are shallow coastal lagoons with diurnal tides and tidal ranges are 0.5-0.7 m (Solis and Powell 1999;Montagna et al. 2012).Long-term estuarine water alkalinity decrease and acidification possibly due to reduced riverine alkalinity export, precipitation decline during drought, freshwater diversion, and calcification in these estuaries, have been reported in a previous study (Hu et al. 2015).To study long-term sulfate variations in estuaries with distinct hydrological conditions, five estuary systems were investigated (Fig. 1; Table 1).One or two major rivers discharging into each estuary were selected based on the methods from Oelsner et al. (2017) and McCutcheon et al. (2019) for the determination of freshwater endmember conditions.

Data availability
Time series estuarine and riverine water sulfate concentration, alkalinity, chlorophyll, pH, and salinity (estimated from chloride concentration) were obtained from Texas Commission on Environmental Quality (TCEQ) website (https://www.tceq.texas.gov/agency/data/wq_data.html) (Table 1).A universal chloride-salinity relationship (salinity (ppt) = 0.00180665 Cl À (mg L À1 )) was used to convert chloride concentration to salinity because of the low detection limit and high dynamic range of the chloride measurement method.Major river discharge data were downloaded from US Geological Survey (USGS, https://maps.waterdata.usgs.gov/mapper/index.html), and the daily mean discharge rates of major rivers, if more than one, that discharge to the same estuary were combined (Table 1).One TCEQ river sampling station, which is less influenced by the tidal cycle, from each major river was selected, and the closest USGS streamgage station was paired (Fig. 1; Table 1).

Conservative mixing model
To eliminate the effect of river-ocean mixing on estuarine water sulfate and alkalinity dynamics, the deviations of sulfate and alkalinity (i.e., ΔSO 4 2À and ΔTA) were calculated by the differences between measured concentrations and theoretical concentrations from their respective mixing models.The mixing models were constructed using river and ocean endmember concentrations.In the sulfate mixing model, river endmember SO 4 2À was directly obtained from the TCEQ river dataset (Table 1) when the corresponding salinity was less than 1, in order to rid the tidal effect.Ocean endmember sulfate was 27451 μmol kg À1 when salinity was 35 (Canfield and Farquhar 2009).In the alkalinity mixing model, river endmember alkalinity was also obtained from the Only one TCEQ estuary station from each estuary system is shown and the complete TCEQ estuary stations can be found in Table S1.TCEQ river dataset when the corresponding salinity was less than 1.Ocean endmember alkalinity was 2427 μmol kg À1 with a salinity of 36.4 (Hu et al. 2015).The slopes and intercepts of each river-ocean mixing line for alkalinity and sulfate were calculated using two points linear regression and then grouped and averaged based on the different river systems discharge to the corresponding estuaries (Fig. 1; Table 1).
River and estuary water quality data after 2000 were included because of the well-documented quality assurance for sulfate measurement (from 1999 to 2000, both EPA 300.0 and SM 4500-SO 4 2À E were acceptable; from 2001 to the present, only EPA 300.0 was acceptable) (Robin Cypher, personal communication).The TCEQ dissolved oxygen data were used to remove sulfate data outliers.When the ratio of ΔSO 4 2À to theoretical sulfate was larger than 10% or less than À10%, the dissolved oxygen level should be over or less than 2 mg L À1 (1.96 mg kg À1 ), respectively, because sulfate enrichment can be only expected under oxic conditions while sulfate removal can only occur under anoxic conditions.In most cases, estuarine water in the study areas was under oxic conditions.The sulfate deviation threshold of 10% was determined from the EPA sulfate measurement uncertainty, and the dissolved oxygen threshold of 2 mg L À1 was from the EPA definition of hypoxia.After outlier removal, only data with ΔSO 4 2À , alkalinity, dissolved oxygen, pH, and river discharge rate available at the same time were further analyzed.In addition, for the chlorophyll-ΔSO 4 2À relationship alone, only data with chlorophyll, ΔSO 4 2À , and dissolved oxygen available at the same time were further analyzed.

Net evaporation estimation and correction
To quantify the effect of net evaporation on the dynamics of alkalinity, freshwater evaporations were first quantified (Appendix SA), and the combined effect of mixing and evaporation on the conservative mixing line was then modeled (Appendix SB).Alkalinity enrichment/removal was quantified using the modified evaporation-mixing relationship.

Salinity normalization approach
Dissolved inorganic carbon (DIC) was calculated from alkalinity, pH, and salinity using CO2SYS.The temperature was set as 22 C; pH scale was set as NBS; Lueker et al. (2000), Dickson (1990), Perez and Fraga (1987), and Lee et al. (2010) were used for carbonic acid dissociation constants, bisulfate dissociation constant, hydrofluoric acid dissociation constant, and borate-to-salinity ratio, respectively.nDIC and nTA were calculated based on the approach in Friis et al. (2003): where TA obs and DIC obs are the observed TA and DIC, respectively; TA 0 and DIC 0 are the zero-salinity TA and DIC determined from linear regression against salinity, respectively; S obs is the observed salinity; S ref is the reference salinity (S ref = 20 in this study).

Results and discussion
Sulfate enrichment in estuaries ΔSO 4 2À scattered around the average theoretical mixing line in all estuaries, and sulfate enrichment (i.e., positive ΔSO 4 2À ) was more pronounced toward the high salinity end in each estuary (Fig. 2).ΔSO 4 2À varied with time, although ubiquitous sulfate enrichment was observed during 2012-2014 (Fig. 3), which coincided with the drought years (Geruo et al. 2017).Besides the very few data points in Sabine Lake, all the other estuaries had little difference in the maximum ΔSO 4 2À and the coresponding time (Fig. 3), even though they had different inflow balances under the normal climate condition (Montagna et al. 2012).The prolonged water residence time during the extreme drought might have minimized the difference in the hydrological conditions and the consequent biogeochemical responses.
The river-ocean mixing line (i.e., the extrapolated river endmember concentration) for sulfate in each estuary did not change much with river discharge variations (Fig. S1), thus applying a universal river-ocean mixing line for each estuary to quantify the nonconservative behavior of sulfate would only cause minimal uncertainties.If there were some river The black line indicates the average theoretical sulfate conservative mixing line calculated using data from all TCEQ river stations.Sulfate enrichment was mostly located toward the high salinity end in each estuary.
endmember sulfate oscillations, the river end sulfate concentrations were still much lower than those commonly found in estuaries (less than 1 vs. over 10 mmol kg À1 ) (Figs. 2, S1).The negligible river endmember sulfate resulted in the evaporation line nearly overlapping with the river-ocean conservative mixing line, and net evaporation would have little effect on the salinity-sulfate relationship compared with the salinity-TA relationship (nonnegligible river endmember alkalinity) (Fig. S1; Appendix SB).The positive ΔSO 4 2À beyond analytical error can only be explained by biogeochemical production.Under drought conditions, the low nutrient input reduced overall autochthonous organic carbon production and allochthonous organic matter input also decreased at the same time.Such conditions could promote autotrophs, especially chemoautotrophs (e.g., sulfide oxidizers) at the sedimentwater interface, to outcompete heterotrophs for oxygen and enrich sulfate.In another estuarine system that experienced drought, Matson and Brinson (1985) found that biochemical oxidation of reduced sulfur in the surface sediment elevated the sulfate/chloride ratio.These authors also found significantly lower chlorophyll concentrations (20% of normal levels), indicating low phytoplankton primary production.However, there was no clear chlorophyll-sulfate relationship in our study areas (Fig. S2), possibly due to the limited time span of chlorophyll data (only before 2006) and the oligotrophic background conditions (N limited) in the south Texas coastal estuaries (Mooney and McClelland 2012).

Potential effect of sulfide oxidation on carbonate chemistry
The redox cycling of sulfur in aquatic environments is closely related to carbonate cycling (e.g., Hu and Cai 2011;Yin et al. 2022;Dias et al. 2023).According to the reaction stoichiometries of several major reactions: 2À across all estuaries except for Sabine Lake, while river discharge rates were low.Color coding indicates the major river discharge rate.Note the differences in colorbar scales.
regardless of the reaction pathway, the expected ratio of ΔSO 4 2À to ΔTA should always be 1 : À2; that is, sulfide oxidation leads to alkalinity loss.During the 2012-2014 drought, the enrichment in sulfate had the potential to decrease alkalinity by up to 35 mmol kg À1 if following this reaction ratio (Fig. 4), and the large alkalinity removal potential was present in all estuaries except for Sabine Lake, where limited data were available.The magnitudes of sulfate enrichment .Expected negative alkalinity deviations were present in all estuaries except for Sabine Lake from 2012 to 2014, and the magnitudes of expected alkalinity deviation were significant compared to observed alkalinity values and variations.Time series alkalinity variations were discussed in Hu et al. (2015).
and depletion make sulfur cycling a significant factor in estuarine water alkalinity dynamics and might contribute to the long-term alkalinity decreasing trend observed from the other study in the same area (Hu et al. 2015).
Estuarine water ΔTA was also calculated based on the estuary-specific river-ocean mixing line, and the relationship between sulfur and alkalinity was further examined using field data (Fig. 5).However, the magnitude of alkalinity reduction was much smaller than expected according to the above reaction stoichiometry, and both positive and negative ΔTA were found even during the extreme drought (Fig. 5A,B).There was a clear latitudinal gradient of ΔTA across all estuaries.For the southern estuaries (Corpus Christi Bay and Aransas Bay), most of the ΔTA were positive, and for the northern estuary (Trinity Bay and Galveston Bay), most of the ΔTA were negative during 2012-2014 (Fig. 5B).The distinct pattern can be caused by variations in mixing model estimation, net evaporation effect, and other biogeochemical processes (e.g., carbonate dissolution) as discussed below.
Despite the fact that reduced sulfur oxidation and alkalinity consumption show close coupling in lab-based incubation experiments (Dias et al. 2023), other factors could also contribute to the observed ΔSO 4 2 -ΔTA relationship in natural conditions.First, the alkalinity levels of the river endmember are often higher than that of Gulf of Mexico seawater, and their variabilities (from 1000 to 7000 μmol kg À1 ) are also higher than river sulfate variabilities (Fig. S1), especially in the southern estuaries.Therefore, a single estuary-specific mixing line for ΔTA calculation would introduce uncertainties, which might lead to the deviation between ΔTA and ΔSO 4 2À (Fig. 5).
Second, groundwater-derived alkalinity could also be an important source of alkalinity in estuaries of the northwestern Gulf of Mexico, especially during drought (Murgulet et al. 2018), although it was not included in the current river-ocean alkalinity mixing model.Groundwater alkalinity fluctuation would complicate the interpretation of alkalinity variations.Multi-endmember mixing models should be considered in the future with concurrent measurement of groundwater-derived alkalinity.
In addition to river endmember variabilities, the observed ΔTA deviations from the conservative mixing line could also be due to the imbalance between evaporation and precipitation when the river endmember value was not zero (Appendices SA, SB).For example, in the southern estuaries within the study area (Aransas Bay and Corpus Christi Bay), evaporation greatly exceeds precipitation (Solis and Powell 1999), and river endmember alkalinity values were much higher than zero (Fig. S1) (see also, Hu et al. 2015;McCutcheon et al. 2019).Net evaporation alone could potentially deviate alkalinity from the conservative mixing line, complicating the interpretation of the observed alkalinity dynamics.In comparison, due to the much lower sulfate concentrations in river endmembers than the oceanic value (Fig. S1), evaporation would have a negligible effect on ΔSO 4 2À .To correct for the evaporation effect on ΔTA, the contribution of net evaporation to freshwater mass balance was first estimated (Appendix SA) and its effect on alkalinity was corrected (Appendix SB).The latitudinal variations in net evaporation affected ΔTA differently in different estuaries (Fig. 5B,C).Nevertheless, net evaporation (up to 30% in Aransas Bay) could only partially explain the discrepancy between ΔSO 4 2À and ΔTA, even in the southernmost estuaries (Fig. 5B,C).Besides, the average net evaporation did not adequately reflect the episodic extreme events.For example, salinity could be over 65 during the extreme drought (Fig. 2), which means net evaporation was at least 44% if the ocean endmember salinity was 36.4,compared to the much smaller net evaporation calculated based on the average residence time (Table A1).To better quantify the nonconservative behavior of alkalinity, more accurate constraints on the conservative mixing model and estimations of net evaporation during episodic extreme conditions are required.However, neither of these issues can be readily addressed using the current dataset.
Estuaries in the northwestern Gulf of Mexico are important habitats for oysters, and a recent study showed that carbonate dissolution and precipitation might occur at the same time along with sulfate enrichment (Dias et al. 2023), and sulfate enrichment (i.e., acid production) could negatively affect oyster fisheries.During the study period, carbonate dissolution, which increases estuarine water alkalinity, and Fe-oxide dissolution, which buffers acidity buildup, may both have contributed to the less than expected alkalinity depletion during the drought (Fig. 5).The relationships between nTA and nDIC (slopes = 0.75-0.98)across different estuaries appeared consistent with the concurrent occurrence of both carbonate dissolution (theoretical slope of 0.5) and sulfide oxidation (no change in nDIC; decrease in nTA) (Fig. S3), although other biogeochemical processes (e.g., C, N, P, Fe, and Mn cycling) cannot be completely ruled out.

Conclusion
Estuaries that experience different hydrological conditions are known for displaying distinct biogeochemical cycling patterns.However, this study revealed common, widespread, sulfate enrichment across all estuaries of the northwestern Gulf of Mexico during an extreme drought from 2012 to 2014.The magnitudes of sulfate enrichment would have caused significant alkalinity removal as sulfide oxidation produces acid.However, such removal was not observed.Sedimentary carbonate dissolution and iron oxide dissolution likely have compensated for this acid production.As the climate becomes drier in the future, the contribution of sulfate enrichment to estuarine alkalinity cycling will likely become more significant, even though the ultimate impact is limited by the inventory of reduced sulfur.To accurately determine the impact of sulfur redox cycling on estuarine alkalinity dynamics over both short and long periods, it is necessary to conduct studies that utilize sophisticated measurements, provide a more detailed characterization of hydrological cycling, and enhance our comprehension of other biogeochemical processes during periods of extreme conditions.

Fig. 1 .
Fig. 1.A map for TCEQ estuary stations, TCEQ river stations, and USGS streamgage stations.The major rivers discharging into each estuary are shown in blue.Only one TCEQ estuary station from each estuary system is shown and the complete TCEQ estuary stations can be found in TableS1.

Fig. 3 .
Fig. 3. Time series of sulfate deviation (ΔSO 4 2À ) across all studied estuaries.From 2012 to 2014, there were distinct positive ΔSO 42À across all estuaries except for Sabine Lake, while river discharge rates were low.Color coding indicates the major river discharge rate.Note the differences in colorbar scales.

Fig. 4 .
Fig. 4. Time series of alkalinity variation and expected alkalinity deviation across all studied estuaries.The expected alkalinity deviation was calculated by two times of -ΔSO 4 2À

Fig. 5 .
Fig. 5. Relationship between ΔSO 4 2À and ΔTA.(A) All data from 2000 to 2020.(B) Data from the drought period (2012-2014).(C) Data from the drought period as in B but with ΔTA corrected for the evaporation effect.Net evaporation corrections were larger in Aransas Bay and Corpus Christi Bay.The red lines in B and C are the theoretical relationship between ΔSO 4 2À and ΔTA considering sulfides and pyrite oxidation and subsequent alkalinity titration.Note the different y-axis scales.

Table 1 .
Identification numbers of TCEQ river and estuary stations, USGS stations, and major rivers discharging into each estuary.*Completestation ID can be found in Supporting Information.