The studied region includes seven estuaries that border the Texas coast in the nwGOM—spanning north to south they are Sabine-Neches Estuary (SNE), Trinity-San Jacinto Estuary (TSJ), Lavaca-Colorado Estuary (LCE), Guadalupe Estuary (GE), Mission-Aransas Estuary (MAE), Nueces Estuary (NE), and Laguna Madre Estuary (LM) (Figure 1). This region hosts 595 km of coastline that lies between 26° and 30°N in a subtropical climate. The nwGOM estuaries are relatively shallow (rarely exceeding 3 m deep), and they are classified as lagoons due to their connectivity and the barrier island system that separates them from the Gulf of Mexico and limits water exchange (Russell and Montagna, 2007; Montagna et al., 2011b). Along with other lagoonal bays that span the majority of the Gulf of Mexico coast, these estuaries make up part of the largest lagoonal system in the world (Dürr et al., 2011).

The area closest to the Louisiana border experiences twice the rate of precipitation as the area closest to the Mexico border. Together with an evaporation gradient, this yields a freshwater balance (sum of all freshwater input via rivers, runoff, and direct precipitation minus evaporation) that is two orders of magnitude higher in the northernmost estuary, i.e., SNE, than in the southernmost estuary, i.e., LM (Montagna et al., 2011a). The differing freshwater balance as well as differences in land use and land cover between estuaries yields differing chemical properties, making these seven estuaries chemically unique despite their physical similarities.

Historical water quality data were obtained from the Texas Commission on Environmental Quality’s (TCEQ’s) Surface Water Quality Monitoring (SWQM) program. TCEQ has collected and compiled water quality data on a quarterly basis across the nwGOM estuaries since 1969, though years of monitoring have varied between stations. We used only those stations with 20 + years of data (including one station in SNE, 10 in TSJ, 13 in LCE, two in GE, five in MAE, 13 in NE, and 10 in LM) and only those stations with 60 + simultaneous measurements for temperature, salinity (or working salinity as defined below), pH, and TA for our analyses. Filtering to remove missing data necessary for carbonate system analyses resulted in removal of all data from 1969 to 1973, the earliest data used in this study were from 1974.

Measurement of TA by TCEQ involved titration to a pH of 4.5 and was monitored using a pH electrode with calibrated acid (H₂SO₄ or HCl) in accordance with EPA Protocol 310.1. This type of interlaboratory methodology would be expected to have a precision of around 5% (Hu et al., 2015). TA is reported by TCEQ in mg-CaCO₃ L^–1, but it was converted to μmol kg^–1 by multiplying by a factor of 20 and dividing by sample density (derived from salinity and room temperature at 22°C) for this analysis. Measurement of pH by TCEQ involved pH sensors that were calibrated daily using a two-point calibration approach with National Institute of Standards and Technology pH 7.00 and 10.00 buffers, and the values were corrected to reflect in situ temperature. Analytical precision was not reported by TCEQ for the pH measurements, although protocol did require measurements to be rounded to the tenth decimal place, so an uncertainty of up to 0.1 can be assumed. Mclaughlin et al. (2017) reported that potentiometric pH sensors in coastal monitoring have error of ± 0.2 pH units, but this is within deployed SeaBird instruments that may have additional sources of error (e.g., much deeper depth). Direct salinity measurements were obtained by TCEQ using a calibrated Hydrolab or YSI sondes using KCl solution of known specific conductance. For those observations where salinity data were missing but specific conductance data were available, a working salinity value was calculated using the linear relationship between specific conductance and salinity (R² = 0.98, as calculated from those records with both salinity and specific conductance data). Working salinity values that were calculated to be ≤ 0 (3% of the 2,747 calculated values) were removed from the analysis.

We used an intensive outlier removal approach that aimed to remove occurrences of falsely reported data due to any human or instrument error coming from the multiple laboratories that analyze TCEQ SWQM samples. A two-step approach for outlier removal based on both bracketed time and salinity intervals was conducted for pH and TA measurements, as was done by Hu et al. (2015). Seven time brackets were established (5-year intervals) and nine salinity brackets were established (with a range of 5 salinity units for salinity < 35 and 35–49.99 and 50 + to encompass the relatively limited observations at hypersalinity). Those observations of each parameter that were deemed outliers within either a time bracket or a salinity bracket were removed. Outliers were defined as measurements that were three times the interquartile range (IQR) above quartile three (Q3) or three times the IQR below Q1 within bracketed subsets of the data. This approach assumes that some outliers are representative of the system and therefore removes fewer data than the standard outlier calculation. Temperature outliers were removed from seasonal subsets of the data using the same outlier designation (eight observations removed). There were no salinity observations that were deemed outliers.

TCEQ’s measured in situ temperature, salinity (or specific conductance that can be converted to working salinity), pH (converted from NBS scale to total scale, see Millero, 2001), and TA were used as inputs to the seacarb package in R for calculations of other carbonate system parameters, propagated error associated with carbonate system parameters, and the β_DIC buffer factor (R Core Team, 2020; Gattuso et al., 2021). For speciation calculations, the carbonic acid dissociation constants K₁ and K₂ were from Millero (2010), the HF dissociation constant was from Perez and Fraga (1987), the HSO₄^– dissociation constant was from Dickson et al. (1990), and boron concentration was from Uppström (1974). Nutrient contributions (phosphate, silica, and ammonium) to TA were not included in speciation calculations, as they were not consistently measured and are assumed to be negligible to TA given the generally low concentrations. Propagated error was calculated using the Gaussian method 1st order, 2nd moment analysis using analytical precision as the standard deviation associated with each of the measured variables (temp = 0.1°C, Sal = 0.1, pH = 0.1, TA = 115 μmol kg^–1).

The calculated Ω_Ar was corrected to account for non-zero calcium concentrations ([Ca²⁺]) in the river end members. Separate corrections were performed for each estuary to reflect different river chemistries. For those estuaries that TCEQ had sufficient [Ca²⁺] data—TSJ, LCE, NE, and LM (59, 26, 58, and 9 aligning [Ca²⁺] and salinity observations, respectively)—a linear regression between TCEQ measured [Ca²⁺] and salinity within each estuary was used to estimate [Ca²⁺] for every date/station observation. A correction to Ω_Ar was then conducted based on the difference between theoretical [Ca²⁺] (with freshwater endmember having zero [Ca²⁺]) and the actual (i.e., estimated) [Ca²⁺], similar to the Ω_Ar correction in McCutcheon et al. (2019). GE, MAE, and SNE all had insufficient [Ca²⁺] data from TCEQ, so Ω_Ar corrections differed. For GE and MAE, similar estuary-wide [Ca²⁺]/Sal regressions were used for correction, using [Ca²⁺] and salinity data from our own lab from 2013 to 2015 and 2014 to 2020, respectively (Hu, unpublished data). For SNE, [Ca²⁺] was estimated for the river endmember by obtaining TA data from TCEQ from the Sabine and Neches rivers (spanning 1974–2013) and dividing TA by two as an estimate of riverine [Ca²⁺] (assuming TA was a result of weathering). The linear relationship between the river endmember (Sal = 0, [Ca²⁺] = 0.28 mmol kg^–1) and the ocean endmember (Sal = 35, [Ca²⁺] = 10.28 mmol kg^–1) was then used to estimate [Ca²⁺], and the data were corrected accordingly. It is worth noting that this Ω_Ar correction method may introduce some level of error (based on fluctuations in endmember Ca²⁺ and Ca²⁺consumption dynamics within each estuary), but the corrected value should be more appropriate than direct seacarb calculations. It has also been demonstrated in NE that Ca²⁺/S relationship can remain linear even when evaporation leads to hypersalinity (McCutcheon et al., 2019).

The buffer factor β_DIC was calculated using the seacarb package Egleston et al. (2010) and the appropriately corrected equation from Álvarez et al. (2014). This buffer factor is particularly relevant since it directly indicates the change in pH that occurs with a given addition of DIC. Given that the common drivers of acidification (e.g., enhanced mineralization and intrusion of atmospheric CO₂) increase DIC, this buffer factor indicates pH response to the acidification mechanism.

Theil-Sen (TS) regression analysis with the Siegel modification (Siegel, 1982) was chosen to generate a series of linear relationships between carbonate system parameters and time (decimal year) at each station; because of its non-parametric approach and robustness to outliers, TS regression is well-suited to explore trends in multi-laboratory agency data (Siegel, 1982; Kaushal et al., 2013; Stets et al., 2014). All analyses were conducted in R, version 4.0.3, and the mblm package was used for regressions (Komsta, 2019; R Core Team, 2020). Rates of increase or decrease in parameters were reported only for those TS regressions that were significant based on alpha = 0.05. Trends are reported at individual stations rather than aggregating into bay units (as was done for investigations of pH and TA trends by Hu et al., 2015) to fully picture spatial heterogeneity. Long-term trends in pH and TA are reported in addition to calculated parameters, as stations were not aggregated into bay entities as was done in Hu et al. (2015) and an additional 5 years of data (2010–2015) were available since the earlier analysis. In addition to pH, we also calculated and reported temporal trends in [H⁺] (Supplementary Table 1); as would be expected, these trends all displayed the inverse sign and had similar levels of significance compared to the respective pH trend (there were only three stations where significance was not found for both trends if found for one, and these stations had marginal significance, i.e., p > 0.01). Based on these regression results and based on the similar average values in pH in the region, we opted to discuss only pH (not [H⁺]), as the pH trends accurately reflected the temporal changes in estuarine acidity (Fassbender et al., 2021).

To confirm that the long-term trends in calculated parameters that were detected by TS regression were valid despite the relatively large propagated error associated with the calculated variables, simulation analyses were conducted to recalculate the slopes of TS regressions over 1,000 simulations at each station for pCO₂, DIC, and Ω_Ar. For each simulation, resampled values for each data point were randomly selected within a normal curve using the calculated propagated error for each calculated parameter as the standard deviation for the curve. A 95% confidence interval of the 1,000 new estimates of the TS regression slope was then calculated; if the confidence interval included a value of zero, then the simulation suggested that there was not a long-term trend in the given parameter. In cases where the simulations showed a significant long-term trend, the mean of all simulated slopes would be the estimate of the long-term trend (see Discussion for the advantages of this approach).

Given the expectation of long-term trends in parameters (and the significant interaction between covariate decimal year and estuaries in analyses of covariance, ANCOVAs), differences between estuaries were investigated using a subset of “current” data. The “current” data consisted of data from (2010 to 2015) for all estuaries except for SNE. Since sampling concluded in 2008 in SNE and sampling frequency was low, data from 1998 to 2008 was used for SNE. Mean differences in carbonate system parameters between estuaries (using current data) were investigated with a single one-way analysis of variance (ANOVA).

In estuaries with sufficient sampling (TSJ, LCE, MAE, NE, LM) stations were categorized as either primary bay stations, i.e., stations in the bay that has direct connection to the GOM, or secondary bay stations, i.e., stations in bays that are separated from the GOM by the primary bays, for spatial comparisons within estuaries (Figure 1). Given that TSJ did not have any sampling in the primary bay areas from 2010 to 2015, the “current” subset was modified to 2000–2005 for only TSJ for primary/secondary comparisons. Differences in parameters between primary and secondary bays were not the same within all estuaries (based on significant interactions in two-way ANOVAs); thus, mean differences between primary and secondary bays were investigated within each estuary using t-tests.

We examined 54 total stations: one station in SNE, 10 in TSJ, 13 in LCE, two in GE, five in, MAE, 13 in NE, and 10 in LM. Salinity showed slight decrease in northern estuaries and mostly increased toward the south (Figure 2A and Supplementary Table 1). Decreasing pH trends dominated in the region (particularly in the mid-coast, NE to TSJ), with rates ranging from −0.0002 to −0.0225 yr^–1 (Figure 2B and Supplementary Table 1). Long-term increases in pH (ranging from 0.0017 to 0.0063 yr^–1) were also observed in one (SNE), one (TSJ), two (NE), and one (LM) stations, respectively (Figures 2B, 4A and Supplementary Table 1). Similarly, decreasing TA dominated in the region (particularly in the mid-coast, NE to LCE), with rates ranging from −1.4 to −30.3 μmol kg^–1 yr^–1. Long-term increases in TA (ranging from 1.9 to 16.7 μmol kg^–1 yr^–1) were also observed at one (TSJ), two (NE), and one (LM) stations, respectively (Figures 2C, 4B and Supplementary Table 1).

Temporal trends in pCO₂ were primarily dominated by long-term increases (particularly in NE, LCE, and TSJ), with nearly half of all stations (25 of 54) demonstrating long-term increases in pCO₂ at rates ranging from 2 to 26 μatm yr^–1; Figures 3A, 4C. Long-term decreases in pCO₂ were also observed (ranging from −2 to −10 μatm yr^–1; Supplementary Table 1), including one (SNE), two (MAE), two (NE), and one (LM) stations, respectively (Figure 3A). Temporal trends in DIC were primarily dominated by long-term decreases (particularly in the mid-coast, NE to LCE), with two thirds of all stations (36 of 54) demonstrating long-term decreases in DIC at rates ranging from −0.9 to −26.1 μmol kg^–1 yr^–1. Long-term increases in DIC were also observed (ranging from 4.7 to 15.0 μmol kg^–1 yr^–1; Figures 3B, 4D and Supplementary Table 1), including one (TSJ), one (LCE), and one (LM) station, respectively.

Temporal trends in Ω_Ar were primarily dominated by long-term decreases (particularly in the mid-coast, NE to TSJ), with nearly two thirds (33 of 54) of all stations experiencing decreases in Ω_Ar at rates ranging −0.0050 to −0.1032 yr^–1 (Figures 3C, 4E and Supplementary Table 1). Long-term increases in Ω_Ar were also observed (ranging from 0.0073 to 0.0292 yr^–1), with increases particularly prominent in LM. Long-term trends in β_DIC were primarily dominated by long-term decreases (particularly in NE, LCE, and TSJ), with more than half (25 of 54) of all stations experiencing decreasing β_DIC at rates ranging from −0.0010 to −0.0050 mmol kg^–1 yr^–1 (Figures 3D, 4F and Supplementary Table 1). Similar to Ω_Ar, long-term increases in β_DIC were prominent in LM; increases were found in seven LM stations as well as one station in each of TSJ, LCE, MAE, and NE, with increases ranging from 0.0010 to 0.0020 mmol kg^–1 yr^–1.

Carbonate chemistry in the region varied substantially between estuaries (Figures 5, 6 and Supplementary Tables 1, 2). SNE, the northernmost estuary, had a significantly different mean for every parameter, including the lowest salinity, pH, TA, DIC, Ω_Ar, and β_DIC, and the highest pCO₂ (Figure 5 and Supplementary Figure 2). Conversely, LM, the southernmost estuary, had the highest salinity, TA, DIC, Ω_Ar, and β_DIC, though only TA and β_DIC were significantly greater than all other estuaries (Figures 5, 6 and Supplementary Figure 2). TSJ also showed relatively unique carbonate chemistry, with relatively low TA, pCO₂, and DIC compared to the rest of the region (Figure 5).

Some of the estuaries showed differences in carbonate chemistry parameters between primary and secondary bay stations (Figure 5). LM was generally the exception to the rule for primary/secondary bay differences across all parameters (Figure 5). With the exception of LM, salinity was lower on average in the secondary bays, although this difference was only significant in LCE (p = 0.0017). pH was lower on average in the secondary bays (only significant in MAE, p = 0.0065). pCO₂ was higher on average in the secondary bay, although this was only significant in MAE (p = 0.0205) and NE (p = 0.0146). Ω_Ar and β_DIC were lower on average in secondary bays, though Ω_Ar was only significantly different in MAE (p = 0.0183) and β_DIC in MAE (p = 0.0478) and TSJ (p = 0.0168). DIC and TA did not vary greatly between primary and secondary bays, although differences were seen in TSJ.

While TCEQ SWQM data are not of the same precision that would be expected in an analytical lab that is specialized in carbonate chemistry studies, the protocols have remained the same over time. Therefore, these data should be adequate for time series analysis, and a number of other studies in the literature have also made this assumption while using government-funded long-term water quality monitoring data (Raymond and Cole, 2003; Raymond et al., 2008; Waldbusser et al., 2011; Stets et al., 2014, 2017; Hu et al., 2015; Barrera and Robbins, 2017; Carstensen et al., 2018; Robbins and Lisle, 2018; Van Dam and Wang, 2019) or even citizen science data (Snyder et al., 2019) for interpreting carbonate chemistry trends.

Propagated error associated with calculated pCO₂, DIC, and Ω_Ar would be expected to be relatively high given the coarse analytical precision of the TCEQ data used in this analysis. Assuming analytical precision of pH as ± 0.1, TA as ± 115 μmol kg^–1 (±5% of mean TA, 2,300 μmol kg^–1), salinity as ± 0.1, and temperature as ± 0.1°C, the mean propagated error associated with calculated pCO₂ values was 184.9 ± 159.1 μatm, the mean propagated error associated with DIC values was 118.7 ± 5.4 μmol kg^–1, and the mean propagated error associated with Ω_Ar was 0.497 ± 0.264.

The simulation method used for validation of long-term trends confirmed the presence and signs of long-term trends in pCO₂, DIC, and Ω_Ar that were found from TS regression with only a few exceptions. The simulations were able to validate the existence and sign of 23 (of 25) of the increasing pCO₂ trends and four (of six) of the decreasing trends. The mean difference between significant slopes calculated by the single TS regression vs. the simulation was −0.2 ± 1.0 μatm yr^–1 (or 3.1 ± 19.3% of the TS slope estimates). The simulations were able to validate the existence and sign of 32 (of 36) of the increasing DIC trends and all three of the decreasing trends. The mean difference between significant slopes calculated by the single TS regression vs. the simulation was −0.3 ± 1.1 μmol kg^–1 yr^–1 (or 4.8 ± 15.1% of the TS slope estimates). The simulations were able to validate the existence and sign of 29 (of 33) of the increasing trends in Ω_Ar and three (of six) of the decreasing trends in Ω_Ar. There was also one station in TSJ (13315) that did not have a significant trend based on the TS method but did have a significant, positive trend based on the simulation. The mean difference between significant slopes calculated by the single TS regression vs. the simulation was −0.0005 ± 0.0051 yr^–1 (or 5.8 ± 21.7% of the TS slope estimates).

While there are some discrepancies between the two methods (single TS regression estimates and 95% confidence intervals of TS regressions of simulated data with incorporated propagated error), the overall picture of the presence and directions of long-term trends in carbonate system parameters within each of the estuaries remains unaltered (Figure 3 and Supplementary Figure 2). This is the first time to our knowledge that the use of calculated carbonate system parameters in the estimation of long-term trends has been validated based on their associated propagated error. It is a promising finding for the utility of government-funded and citizen science-sponsored data, and it provides support for the continued sponsorship of such programs for understanding changes in coastal carbonate chemistry.

It is worth noting that DIC and pH are among the better combinations of parameters for carbonate system calculations (Hunt et al., 2011; Patsavas et al., 2015). Since the TCEQ data did not include measures of DIC, TA and pH were used for speciation calculations; this simple use of TA has been shown to lead to overestimation of pCO₂ (Hunt et al., 2011) since organic alkalinity is not accounted for. Organic alkalinity has been assumed to be negligible in other studies when only TA and pH were available for speciation calculations (Prasad et al., 2013; Carstensen et al., 2018). Present day organic alkalinity in the region (MAE and NE, calculated by the difference in measured TA and TA calculated with DIC and pH) usually does not exceed 50 μmol kg^–1, or less than 2% of TA (Yao and Hu, 2017; McCutcheon et al., 2019). It is assumed that there are not long-term trends in organic alkalinity. However, analysis based on the agency collected data suggest that total organic carbon concentrations in many nwGOM estuaries likely exhibit long-term decline (Montagna and Palmer, 2012; CTBS, 2017). It is reasonable to assume that organic alkalinity level would show a positive correlation with total organic carbon, hence the pCO₂ trends calculated based on the assumption of no organic alkalinity trend would represent the lower estimates of “true” values.

Around the world, several estuaries have now been shown to have long-term trends in pH and TA (Table 1). with a meta-analysis revealing equal likelihood of increasing or decreasing pH trends (Carstensen and Duarte, 2019). Baumann and Smith (2018) investigated pH trends across 16 U.S. Nation Estuary Research Reserve System sites. While interannual pH change rates are not provided, they reported that seven of the 16 sites have decreasing pH, four have increasing pH, and the magnitude of change is ∼10 times that of the open ocean (∼−0.0017 yr^–1, González-Dávila et al., 2007; Byrne et al., 2010). The magnitude of annual change in pH and TA that we have reported is mostly in line with trends that have been reported in other estuaries (Table 1). Only one other study that we are aware of has reported long-term trends in calculated carbonate system parameters; in Danish Fjords, Carstensen et al. (2018) reported mostly long-term increases in pCO₂ and DIC (Table 1). While such increasing DIC would generally be expected in oceanic environments with increasing pCO₂ (both resulting from CO₂ uptake), the widespread decreasing trends that were found here likely directly result from the long-term decreases in TA in the nwGOM region. Carstensen et al. (2018) did not report values of annual change in Ω (although they did note increasing Ω and decreasing Ω surrounding a regime shift in one estuary). The large spatial variation in carbonate chemistry within estuaries is often due to control exerted by NEM. Heterogeneity in NEM becomes apparent at distances of less than 2 km in the nwGOM estuaries (Russell and Montagna, 2007). More freshwater-influenced regions of estuaries are generally more net heterotrophic than the rest of the estuary due to allochthonous carbon input (Caffrey, 2004) and river water itself is generally supersaturated with CO₂ (Butman and Raymond, 2011; Joesoef et al., 2015); thus, the elevated pCO₂ and depressed pH, Ω_Ar, and β_DIC found in the secondary bays relative to the primary bays (Figure 5) are spatial relationships that could be noted along the freshwater gradient within many estuaries.

Like many other coastal areas (e.g., Scanes et al., 2020; Da et al., 2021), the nwGOM estuaries have experienced gradual but significant warming over the past decades, with the overall warming rate of 0.0428°C yr^–1 (Montagna et al., 2011b). This change alone would lead to 1–2 μatm yr^–1 increase in pCO₂ and ∼0.0004–0.0006 yr^–1 decrease in pH, based on the average water chemistry parameters and an average temperature of 25°C in Supplementary Table 2. It appears that warming alone would explain a small fraction of the observed trends, as carbonate equilibrium would shift toward higher pCO₂ and lower pH (Table 1). However, warming can also lead to enhanced aquatic respiration (e.g., Hopkinson and Smith, 2005; Ni et al., 2019), which can further increase pCO₂ and tilt the NEM balance toward heterotrophy. There is evidence of this shift in NEM in a portion of NE (Corpus Christi Bay, based on a 1982–2002 dataset), where the decrease in surface dissolved oxygen (DO) concentration exceeds that which can be explained by thermally induced oxygen solubility decrease (Applebaum et al., 2005). Similarly, Galveston Bay also exhibits declining DO in its subbays (Gonzalez, 2011) as this highly urbanized estuary undergoes eutrophication (Bugica et al., 2020).

For the few sites that showed decreasing pCO₂ over time, the SNE, NE, and LM sites all had increases in pH and decreasing TA (or DIC for the LM site), while the two MAE sites with decreasing pCO₂ trends both had decreasing TA and either decreasing or no trend in pH (Supplementary Table 1). A change in hydrological conditions might explain these trends, as reduction of alkaline freshwater leads to TA decrease and presumably freshwater has higher pCO₂ than the ocean water (Hu et al., 2015). Meanwhile, reduced freshwater input could hamper the delivery of organic matter, hence the pCO₂ trends at these stations indeed reflected the lower estimates because of the assumption of non-variant organic alkalinity in the speciation calculations. Our ongoing study at MAE may shed light on this issue in the future with better quality carbonate chemistry characterization.

Trends in pCO₂ have implications for air-water CO₂ flux (Landschützer et al., 2013). Given that air-water CO₂ flux is a function of the gradient between atmospheric and surface water pCO₂, long-term trends in pCO₂ of surface water that differ from long-term trends of atmospheric CO₂ indicate changes in the water’s magnitude of a carbon sink (absorbs atmospheric CO₂) or carbon source (outgasses CO₂ to the atmosphere) (Landschützer et al., 2013). For comparison to pCO₂ trends in estuaries, the average annual change in atmospheric pCO₂ between 1974 and 2015 was 1.72 μatm yr^––1.¹ Given that the majority of the stations throughout the nwGOM estuaries have experienced increasing pCO₂ exceeding the atmospheric rate (Figures 6C,D, and Supplementary Table 1), and given that estuaries tended to have pCO₂ higher than that of the atmosphere (Supplementary Table 1), we can safely conclude that the nwGOM estuaries have increased in their level of CO₂ outgassing over the observed decades.

Rates of increase in pCO₂ were generally higher across all nwGOM estuaries than surface ocean pCO₂ trends (1.88 μatm yr^–1; Doney, 2010). This elevated pCO₂ increase relative to the global ocean has also been observed across the entire nwGOM coastal ocean (3.20 μatm yr^–1; Kealoha et al., 2020). Individual estuarine stations that had increasing pCO₂ in this study nearly always exceeded the rate found for the nwGOM coastal ocean, which would indicate that advection of estuarine waters into the coastal GOM may increase coastal pCO₂ trends. However, given the lack of significant long-term trends across all estuarine stations, it is possible that better estuarine sampling coverage would reveal a spatially weighted mean similar to or less than the adjacent coastal ocean.

CO₂ flux estimates from recent years in the nwGOM estuaries have been shown to vary between estuaries, with MAE and NE acting as mild CO₂ sources (6.9 ± 6.5 and 2.7 ± 8.1°C m^–2 d^–1, respectively) and GE and LCE acting as moderate CO₂ sources (35.9 ± 24.2 and 25.6 ± 28.8°C m^–2 d^–1, respectively) on an annual basis (Yao et al., 2020). Wet/dry cycles play a substantial role in the interannual variability of CO₂ fluxes in the region (Yao et al., 2020). As pCO₂ trends have varied between increasing and decreasing in NE and MAE, the flux in there may not have changed drastically over the past decades. However, pCO₂ is clearly increasing in LCE, indicating that it likely had lower CO₂ fluxes decades ago. Yao et al. (2020) did not consider the two northernmost (SNE and TSJ) nor the southernmost (LM) estuaries. Due to SNE and LM experiencing the extremes of substantial precipitation and evaporation gradients, their carbonate system characteristics differed greatly from the other estuaries in the nwGOM (Figure 5), as did the long-term trends that they have experienced (Figure 4 and Supplementary Table 1). It is important to understand long-term trends such as these when fluxes are applied to regional or global carbon budgets.

Trends in Ω_Ar have implications for biological calcification as well as calcium carbonate preservation and burial (Fabry et al., 2008; Waldbusser et al., 2015). The nwGOM estuaries is inhabited by Eastern Oysters, ecologically and economically important species that require suitable carbonate chemistry to build and maintain their shells (Buzan et al., 2009; Pollack et al., 2011; Beseres Pollack et al., 2013). With the exception of SNE, nwGOM estuaries have generally remained supersaturated with respect to aragonite (Supplementary Table 2), but even slight decreases in Ω_Ar while still maintaining supersaturation have been shown to negatively affect calcifying organisms (Gattuso, 1998). Trends in Ω_Ar in the nwGOM estuaries were generally decreasing, except for LM (Figure 4). Given that Ω_Ar still remains well above saturation (except for SNE), carbonate chemistry is likely a secondary stressor to changes in salinity in the region. Oyster populations have always been more abundant in the lower-salinity estuaries of the Texas coast (Montagna et al., 2007), but salinity has been increasing along much of the Texas coast over past decades (Bugica et al., 2020), potentially threatening oysters with increased predation, increased disease, decreased settlement, and salinity tolerance that is further compromised by increasing water temperature (VanderKooy, 2012). Despite the more pressing threat, this declining trend in Ω_Ar may have already reduced biological calcification rates over past decades, and if the trend continues it may threaten the shell budget and subsequently threaten sustained oyster populations into the future (Powell and Klinck, 2007. It is also worth noting that the estuary that has experienced the greatest rate of decrease in Ω_Ar (TSJ, Figure 4E) is also the estuary with the heaviest reliance on commercial oyster harvests.

Trends in buffer factors have implications for fluctuations in other carbonate system parameters, and β_DIC is particularly indicative of changes in pH. As noted by Hu et al. (2015), those stations with decreasing pH had a much higher rate of decrease than that in the open ocean, which is around 0.0017 yr^–1 (González-Dávila et al., 2007; Byrne et al., 2010; Supplementary Table 1). The long-term decreases that have been observed in β_DIC in the nwGOM estuaries over the past several decades indicate that the pH response of these systems has become increasingly sensitive to DIC additions. Therefore, short-term processes like the stratification and development of hypoxia that can occur in the region (Applebaum et al., 2005; McCutcheon et al., 2019) or long-term processes like decreasing freshwater inflow or increasing atmospheric CO₂ all have the potential to cause greater fluctuations and disruption to the pH now than the same processes would have decades ago. This long-term decrease in buffer capacity and increase in potential for change has been described in ocean environments previously (Shaw et al., 2013), but it has never explicitly been noted in estuarine environments. Despite this decrease in β_DIC, β_DIC is generally high in nwGOM estuaries, often exceeding the values in nwGOM open water, ∼0.26–0.29 mmol kg^–1 (Figures 6E,F) when buffer capacity is generally lower in estuarine waters.