LIS is an urbanized, seasonally eutrophic estuary on the east coast of the United States ( Figure 1 ). The macrotidal estuary is bordered by Connecticut to the north, by Long Island, NY to the south, and by New York City to the west. Most of the seawater exchange occurs through the mouth of the estuary in eastern LIS (ELIS) where it is connected to the Mid Atlantic Bight (O’Donnell et al., 2014; Whitney et al., 2016). A much smaller exchange also occurs in WLIS where the head of the estuary is connected to the East River, a tidal strait connected to the Hudson River estuary (Gay et al., 2004). The Connecticut River is the largest freshwater source, accounting for 75% of annual freshwater delivery, and enters near the estuary mouth. The next largest freshwater source is the Housatonic River in western LIS, accounting for 10% of annual river discharge to the estuary (Gay et al., 2004). Due to the circulation patterns in LIS, the Connecticut River comprises >50% and >70% of the freshwater in western and eastern LIS, respectively, while the Housatonic River comprises >25% and >10% of the freshwater in western and eastern LIS, respectively (Deignan-Schmidt and Whitney, 2018). The residence time of freshwater in LIS ranges from 3-6 months, depending on river discharge (Deignan-Schmidt and Whitney, 2018).

Carbonate system, nutrient, and hydrographic parameters were observed at ten stations across thirteen sampling events between March 2020 and August 2022 aboard the R/V John Dempsey ( Supplementary Table S1 ). Each sampling event spanned at least three days, with one day each dedicated to the eastern, western, and central Sound, normally in that order. The sampling frequency is always greatest in the summer (typically biweekly) to allow CT DEEP to better assess patterns of hypoxia, with monthly sampling occurring year-round. Thus, a full seasonal comparison between spring, summer, and autumn was captured for 2020 and 2021.

Discrete samples for TA and DIC were collected at each station in Figure 1 for both surface (2 m depth) and deep waters (~2 m above bottom). DIC samples were collected in triplicate from 10L Niskin bottles into 20 mL borosilicate vials with polytetrafluoroethylene (PTFE) and aluminum crimp tops and preserved with 10μL saturated mercuric chloride to a final concentration of 0.05% (v/v). DIC was measured using a Shimadzu TOC-L analyzer (USGS, 2019) calibrated using CO₂ certified reference materials (CRMs) (Dickson et al., 2003) with an average precision of 4.76 μmol kg^-1 (< 0.3%). TA samples were collected in 250 mL high density polyethylene (HDPE) bottles and preserved with 100 μL of saturated mercuric chloride. TA was analyzed using a Contros HYDROFIA TA Analyzer (Seelmann et al., 2019) calibrated with CO₂ CRMs with an average precision of 2.5 μmol kg^-1 (< 0.2%). Orthophosphate and silicate samples were collected and analyzed by CT DEEP following EPA methods 365.1 and 370.1, respectively (EPA, 1978, 1993).

Temperature (± 0.01°C), salinity (practical salinity scale, ± 0.1), dissolved oxygen (DO) (± 3.1 μmol kg^-1), pH (NBS scale, ± 0.1), and partial pressure of CO₂ (pCO₂) (± 2-5 ppm) were measured across LIS surface waters. These parameters were measured continuously in seawater pumped onto the ship from a ~2-m depth intake. Temperature, salinity, and pH were measured once per minute using a YSI Exo2 multiparameter sonde. A three-point pH calibration (4.01, 7.00, and 10.01, all ± 0.01) and one-point salinity calibration (50 mS cm^-1 conductivity, Ricca) were performed at 25°C before each cruise. pCO₂ was measured once per second in the flow-through surface water using a Turner C-Sense pCO₂ probe and then averaged across 5-minute intervals. The probe was calibrated annually by the manufacturer and displayed no drift over the study period.

Carbonate system computations including error propagations were carried out in Python using the PyCO2SYS package (Humphreys et al., 2022) with inputs of TA and DIC along with salinity, temperature, and pressure. The empirical constants used were those of Lueker et al. (2000), Dickson (1990), Dickson and Riley (1979), and Lee et al. (2010). For all cruises labelled “WQ” (water quality, Supplementary Table S1 ), nutrient data were available from CT DEEP and phosphate and silicate were included in calculations. However, nutrient observations were not available for “HY” (hypoxia, Supplementary Table S1 ) cruises, so concentrations were assumed to be 0. This assumption was made based on high rates of surface nutrient uptake in the summer, as well as the relatively low error associated with omitting the alkalinity contribution from nutrients (0.5 to 1.5 μmol kg^-1), which was propagated into the calculations of the carbonate system.

Apparent oxygen utilization (AOU) was calculated in Python using the Gibbs Seawater package (McDougall and Barker, 2011). In-situ salinity, temperature, and depth data were used to calculate potential temperature, from which the oxygen solubility was calculated via the equations of Garcia and Gordon (1993).

The precipitation conditions in Connecticut, and more broadly the northeast US, are highly relevant to the hydrological dynamics of LIS, as > 85% of freshwater input to the estuary occurs from rivers traversing the state. Figure 2 shows the 2020 – 2022 monthly average discharge of the two major tributaries to LIS, the Housatonic River in the center/west (accounting for 10% of total freshwater input to LIS), and the Connecticut River in the east (accounting for 75% of total freshwater input to LIS, Gay et al., 2004). The monthly river discharge climatology shows that in summer 2021, the discharge of both rivers exceeded that of 2020 and 2022 by 3-4 times. This was reflected in the salinity of WLIS, as summertime (August) salinity was ~1 unit lower in 2021 than in 2020 and 2022 ( Supplementary Figure 1 ).

This study captured the contrast between two historic drought seasons (summer and fall of 2020 and 2022) and one high-precipitation season (summer 2021). In August and September 2020, up to 60% of the state of Connecticut experienced a standardized precipitation index (SPI) < -1.6, indicating extreme drought conditions relative to the climatological mean for the region (NIDIS, 2023). SPI remained < -0.5 for much of the state from July 2020 through April 2021. During May 2021 – February 2022, exceptionally wet conditions prevailed, with ~40% of Connecticut experiencing an SPI > 2 (NIDIS, 2023). Four tropical storm systems made landfall in Connecticut between July and September of 2021 (Elsa, Fred, Henri, and Ida), accounting for much of the excess rainfall. The wet conditions of 2021 were followed by an even more extreme drought in 2022, where from April to August 60% of the state had an SPI < -1.6, and 20% < -2.0 (NIDIS, 2023).

For all sampling events, TA was generally highest in ELIS at the estuary mouth and decreased westward to the estuary head ( Supplementary Figure 2 ), consistent with salinity trends (p < 10^-35, R² = 0.59, n = 178, Figure 3A ) and seawater exchange. Throughout the study period, surface and bottom water TA ranged from 1862 – 2113 and 1840 – 2133 μmol kg^-1, respectively. For all three sampling years, the highest TA values tended to occur in the fall (late August – October). Summer and autumn TA values were lower in 2021 as compared to 2020 and 2022, especially in surface waters. The maximum surface TA value in August/September 2021 was 2017 μmol kg^-1, as compared to 2070 and 2069 μmol kg^-1 in 2020 and 2022, respectively. In eastern LIS (ELIS), the influence of increased Connecticut River discharge in September 2021 is clear, as the region near the river mouth was locally depleted in TA, with a value of ~1875 μmol kg^-1 as compared to surrounding values 50-100 μmol kg^-1 higher.

Because DIC includes the concentration of CO₂, it is more strongly affected by biological activity than TA. WLIS is highly productive in the summer, so DIC displayed more spatial and seasonal variation than TA, but still correlated strongly with salinity (p < 10^-21, R² = 0.38, n = 200, Figure 3B ). During the study period, DIC in surface and bottom waters ranged from 1678 – 1908 and 1704 -1984 μmol kg^-1, respectively, generally increasing with depth and towards the estuary mouth. The vertical gradient in DIC had a strong seasonal variation, with the largest gradient occurring in the stratified summer months and the smallest gradient in the fall ( Supplementary Figure 3 ). In the summer, high productivity in the surface waters depressed DIC, while comparatively, DIC was elevated in deeper waters, consistent with intense stratification and subsurface respiration of organic matter. As with TA, the interannual variation of DIC was dominated by autumn minima in 2021 as compared to 2020 and 2022. In October 2020 and August 2022, surface DIC values were ~1850-1900 μmol kg^-1. By contrast, in September and October of 2021, DIC was always less than 1800 μmol kg^-1 and was as low as 1680 μmol kg^-1 in WLIS.

Aragonite saturation state (Ω_ar) was calculated from TA and DIC measurements, yielding a mean propagated error of 0.04 units. Ω_ar is an important acidification metric to consider for calcifying organisms sensitive to carbonate ion concentrations (Waldbusser et al., 2015), though species-specific thresholds are not well-defined in this region. In general, Ω_ar < 1 (aragonite undersaturation) is considered deleterious to calcifying organisms, as aragonite dissolution is thermodynamically favored. However, Ω_ar < 1.5 or higher can still be a detrimental condition depending on the organism and its life stage (Ekstrom et al., 2015). The trend in Ω_ar ( Supplementary Figure 4 ) closely followed that of DIC, with a strong summertime along-estuary gradient in surface waters. Ω_ar ranged from 0.87 – 3.01 and 0.59 – 2.77 in surface and bottom waters, respectively. Notably, the lowest Ω_ar values were observed during summer 2020, during which much of WLIS bottom waters were aragonite undersaturated for at least three sampling events spanning six weeks in July and August. Ω_ar was low, but exceeded 1, in WLIS during summer 2022. In summer 2021, a year with anomalously high precipitation rates, bottom water Ω_ar always exceeded 1.3. Despite lower summertime estuary-wide TA due to increased freshwater discharge, Ω_ar remained oversaturated in bottom waters due to concurrent low DIC, and thus a low ratio of DIC to TA.

Assuming that Ω_ar is primarily biologically mediated (via photosynthesis/respiration and calcification/dissolution), we expect a strong relationship with apparent oxygen utilization (AOU), especially during the stratified summer season. AOU describes the difference between the observed dissolved oxygen (DO) concentration and that at saturation. Positive and negative AOU values indicate DO under- and over-saturation, respectively. The use of AOU as a proxy for photosynthesis and respiration has the limitation that oxygen equilibrates rapidly in the surface mixed layer and can be advected from adjacent water parcels, so the AOU of a given water parcel does not solely reflect oxygen production or utilization in that location. However, AOU is still useful for a coarse estimate of biological oxygen consumption and production. Figure 4 shows the relationship between Ω_ar and AOU during this study period. While the two parameters showed strong correlation in 2020 (R² = 0.70, p < 10^-16, n = 68) and 2022 (R² = 0.51, p < 0.001, n = 23), the relationship was less significant during 2021 (R^{2 =} 0.07, p < 0.10, n = 48). Longer estuarine residence times enhance subsurface oxygen depletion and CO₂ accumulation (e.g., Vlahos and Whitney, 2017; Zeldis et al., 2022), which likely explains the reduced correlation between AOU and Ω_ar during high-discharge conditions leading to shorter residence time. For the entire dataset, the relationship between Ω_ar and AOU was significant (R² = 0.26, p < 10^-9, n = 139).

The ratio of DIC to TA (DIC/TA) is another proxy for acidification that is more straightforward and precise than Ω_ar, as it relies on simple division of two measured parameters rather than calculations using empirical constants. DIC consists of CO₂, bicarbonate, and carbonate, while TA consists of carbonate (contributing twice as much to TA as to DIC), bicarbonate, and other proton acceptors. Thus, higher DIC/TA indicates a greater proportion of CO₂ and thus more acidified waters, while lower DIC/TA indicates a more buffered water mass. For context, the seawater endmember of DIC/TA is ~0.87 (Cai et al., 2020). In this study, DIC/TA ranged from 0.86 (WLIS surface, September 2021) to 1.00 (ELIS bottom, July 2020). The highest DIC/TA values (most acidified) were generally observed during 2020 (0.88 – 1.00, median 0.94, n = 69) and the lowest values (least acidified) during 2021 (0.86 – 0.95, median 0.91, n = 48), with 2022 falling in the middle (0.90 – 0.96, median 0.92, n = 23), indicating that LIS was overall better buffered during the high-precipitation summer (2021) than in the drought summers (2020 and 2022). DIC/TA and Ω_ar were strongly negatively correlated for the entire study period ( Figure 4B , R² = 0.98, p < 10^-100, n = 139).

Hypoxia (defined as DO < 3 mg L^-1 or ~94 μmol kg^-1 in LIS) has decreased in duration and extent due to mandated reductions in nitrogen emissions to the estuary (Whitney and Vlahos, 2021; Long Island Sound Study, 2022). The hypoxic duration and extent in 2020 (43 d and 163 km², respectively) were among the lowest observed in LIS since the implementation of the nitrogen TMDL in 2000. In 2021, hypoxic duration and extent were 47 d and 368 km², respectively, while in 2022, hypoxic duration and extent were 57 d and 225 km², respectively. The hypoxic area of LIS during all three years represented in this study was below the 539 km² average for the entire time series from 1991-2021. Table 1 summarizes these trends.

In many cases, estuarine eutrophication leads to concurrent low-DO and high-CO₂ conditions (Cai et al., 2011; Wallace et al., 2014; Shen et al., 2022), though these conditions may be temporally and spatially decoupled (Shen et al., 2023). While the most intense hypoxia was consistently observed at WLIS station A4 in August for the three study years, the co-occurring acidification trends were more variable. In 2020, a drought summer, though TA and thus Ω_ar data were not available for bottom waters during August, the lowest DO observed (55.3 μmol L^-1, station A4, WQAug20) was concurrent with a bottom-water DIC maximum (1984.2 μmol kg^-1) and pH_nbs minimum (7.6). In 2021, the high-precipitation summer, the lowest observed DO (58.1 μmol L^-1, station A4, HYAug21) was similar to the concentration observed in 2020, but DIC was much lower (1839.7 μmol kg^-1) while Ω_ar was 1.79, well above saturation conditions. Finally, in 2022, the lowest observed DO (44.7 μmol L^-1, station A4, WQAug22) occurred with low Ω_ar (1.09) and elevated DIC (1946.2 μmol kg^-1). Supplementary Figure 5 shows the cross-sectional distribution of DO for each sampling event. In summary, the co-occurrence of hypoxia and acidification was clear for times of drought in 2020 and 2022, but not for the summer with high precipitation (2021).

Though the spatial extent of this study captured both western and eastern LIS, this discussion primarily focuses on western LIS (WLIS), as that is the region where the most severe hypoxia and acidification have been observed.

Of the three study years, the carbonate system dynamics of 2020 and 2022 were comparable, while 2021 showed a divergent pattern. 2020 and 2022 were both drought years, with especially low freshwater discharge during the summer, while in 2021 there was increased freshwater input, primarily due to tropical storms. In some estuaries, lower river discharge is associated with reduced hypoxia, as reduced freshwater flow leads to less nutrient loading, lower primary productivity, and less organic material fueling microbial respiration (Zheng et al., 2016; Shen et al., 2022). By contrast, in LIS low-flow years have been shown to be net heterotrophic, as longer estuarine residence time leads to DO depletion, despite the reduced source of allochthonous organic material (Vlahos and Whitney, 2017). Despite the slightly greater hypoxic area observed in 2021 (368 km²) relative to 2020 and 2022 (163 and 225 km², respectively), the entire study period had hypoxic conditions well below the mean area (538 km²) for the monitoring period 1991-2021 (see Table 1 ). However, the westernmost sampling station A4 still displayed hypoxic conditions in bottom waters for several weeks of monitoring during each sampling year. Thus, freshwater discharge conditions did not appear to substantially impact summertime hypoxia in LIS during the study period.

During the timing of minimum DO concentration (July-August), hypoxic conditions were relatively similar in the bottom waters of WLIS ( Supplementary Figure 5 ) for all three years, but DIC was nearly 100 μmol kg^-1 higher in the drought summers 2020 and 2022 as compared to the wet summer in 2021 ( Table 1 ; Supplementary Figure 3 ). Ω_ar was also substantially higher in 2021 than in 2020 and 2022 ( Table 1 ; Supplementary Figure 4 ). This indicates that the variation in freshwater discharge had more influence on bottom-water acidification than on hypoxia. We propose that these differences are due to differences in the chemistry of the freshwater endmember, variable proportions of fresh- and seawater endmembers, the differing equilibration times for oxygen and CO₂, and variability in estuarine residence time.

DIC and TA data are available for freshwater endmembers in the Connecticut and Housatonic Rivers during 2020 and 2021 (McGuinness, 2022), overlapping with two of the three years in this study. These two rivers account on average for ~55 and 25% of the freshwater content in WLIS, respectively (Deignan-Schmidt and Whitney, 2018). The Housatonic River, which discharges near the eastern boundary of WLIS, has relatively higher TA (averaging 1578 ± 236 μmol L^-1 during 2020-2021) than the Connecticut River, likely due to the presence of carbonate rock formations in the watershed (Robinson and Kapo, 2003). However, the main geologic formation of the Housatonic watershed is composed of slates and pelitic metasedimentary rocks. By contrast, the Connecticut River has relatively lower TA (averaging 531 ± 110 μmol L^-1 during 2020-2021) as it primarily flows through granite and gneiss bedrocks, though there are some carbonate-containing formations in central Connecticut (Robinson and Kapo, 2003). WLIS may have a more substantial source of freshwater alkalinity, especially under high-flow conditions, as it contains a greater proportion of Housatonic River outflow than ELIS.

While river TA was variable among sampling events, there was no statistically significant relationship with discharge (Connecticut: R² = 0.06, p = 0.52, n = 8, Housatonic: R² = 0.05, p = 0.59, n = 8). However, river DIC was negatively correlated with discharge for 2020-2021 (Connecticut: R² = 0.72, p = 0.004, n = 9, Housatonic: R² = 0.64, p = 0.02, n = 8), and thus DIC/TA was also negatively correlated with discharge (Connecticut: R² = 0.45, p = 0.07, n = 8, Housatonic: R² = 0.81, p = 0.01, n = 6). River DIC likely decreased under elevated river discharge because of a shorter residence time for organic matter respiration and CO₂ accumulation.

The increased buffering capacity of the Housatonic and Connecticut Rivers under high discharge conditions, as indicated by decreased DIC/TA, influenced the acid-base chemistry of WLIS during summer/fall 2021 ( Figure 5 ). Though bottom-water TA was relatively consistent in WLIS during the study period, DIC and thus DIC/TA were much lower in 2021 than in 2020 and 2022. This explains the relatively higher Ω_ar (> 1.5) in WLIS bottom waters during peak hypoxia in 2021, relative to the undersaturated (Ω_ar < 1) and near saturated (Ω_ar < 1.1) conditions in the drought summers of 2020 and 2022, respectively. In 2012, Wallace et al. (2014) also showed summertime aragonite undersaturation in the bottom waters of WLIS co-occurring with hypoxic conditions. Hypoxic extent was relatively high (747 km²) in 2012 relative to the 1991-2021 average of 538 km². 2012 was also a historic drought year in Connecticut, on par with the conditions of 2020 and 2022. Though river chemistry data is not available for the 2012 study, it is likely that low river discharge increased the DIC/TA value of the freshwater endmember, reducing the buffering capacity of WLIS. The increased water residence time associated with drought likely also contributed to the depletion of O₂ and the accumulation of CO₂. Our observations during a high-flow year (2021) indicate that the conditions leading to hypoxic conditions in LIS do not always simultaneously induce acidification. The buffering capacity of the incoming freshwater appears to play a role in mitigating or exacerbating summertime acidification.

Though precipitation in the northeast US is expected to increase under climate change, most of this increase is predicted to occur in the winter and spring, with little change in the summer (USGCRP, 2018). By contrast, late summer and fall streamflow are expected to decrease due to a greater chance for summertime drought (Demaria et al., 2016). If these regional projections hold true for the LIS watershed, a decrease in alkalinity loading from the Housatonic River and an increase in freshwater DIC may exacerbate late summer acidification in WLIS.

The unique hydrodynamics of LIS can complicate comparisons to other estuaries. While most US east coast estuaries have deltaic features, with the main freshwater source at the head mixing gradually with seawater towards the mouth, LIS was formed through the glacial deposition of Long Island, and the largest freshwater source (Connecticut River) is near its mouth. Despite these unique physical characteristics, valuable comparisons can be made regarding the summertime stratification, hypoxia, and acidification displayed in other major US estuaries including the Delaware Bay and Gulf of Mexico. These estuaries are discussed due to their data availability and comparable dynamics to LIS.

Delaware Bay is a drowned river valley or coastal plain estuary which is a tributary of the mid-Atlantic Bight and receives its largest freshwater input from the Delaware River at the estuary head. Delaware Bay has a similar seasonal weather pattern as LIS, and the chemistry of both the fresh- and seawater endmembers are relatively similar, though TA is higher in the Housatonic River (1578 ± 236 μmol kg^-1) than in the Delaware River (~750-1000 μmol kg^-1, Joesoef et al., 2017). An analysis of the relationship between seasonal discharge and the resulting acidification of Delaware Bay showed that high discharge decreased both DIC and TA and increased the DIC/TA ratio (Joesoef et al., 2017). This directly contrasts with observations in LIS, wherein the DIC/TA ratio was lower under high-discharge conditions due to a decrease in DIC, but no significant change in TA. Joesoef et al. (2017) suggested that DIC/TA increased under high-flow conditions in Delaware Bay because of enhanced uptake of soil organic matter and reduced CO₂ outgassing. This contrast points to differing hydrodynamic and biogeochemical processes in LIS than those observed for Delaware Bay, which should be investigated in a pan-watershed approach to further understanding of LIS acidification dynamics.

The Louisiana Shelf of the Gulf of Mexico undergoes seasonal hypoxia due to agricultural nutrient pollution in the Mississippi River watershed (Rabalais and Turner, 2019). While subsurface acidification is prevalent in these coastal waters (Wang et al., 2020), the high alkalinity of the Mississippi River (~2400 μmol kg^-1, Savoie et al., 2022) provides resilience against acidification. Despite similar eutrophic and low-oxygen conditions in the Gulf of Mexico as in LIS, the associated changes to pH and Ω_ar are less dramatic than in LIS due to alkaline river water input and warmer waters that hold less CO₂ (Cai et al., 2020). Due to the high alkalinity of the Mississippi River, Savoie et al. (2022) observed that under high-flow conditions (opening of a spillway), surface water TA and Ω_ar increased in the Mississippi Sound, while subsurface Ω_ar decreased in association with stratification and hypoxia. In our observations in WLIS associated with increased Housatonic River discharge, we observed similar surface water dynamics. However, TA and Ω_ar were higher throughout the water column in WLIS. This divergence may be attributed to the greater depth of LIS (~40m compared to 4m in the Mississippi Sound) and its partially well-mixed characteristics, especially during the multiple storm events of summer 2021 which brought both increased precipitation and winds.

Comparisons of LIS to Delaware Bay and the Gulf of Mexico provide context to inform future studies, despite clear differences between the estuaries. For example, the decrease in bottom water DIC/TA under conditions of high freshwater input appears unique to western LIS among these estuaries.

In this study, we demonstrated seasonal and interannual variations in carbonate chemistry dynamics through observations of TA and DIC across LIS from 2020-2022 on monthly to seasonal intervals. Both TA and DIC were highly correlated to salinity and increased from the head to the mouth of LIS. Seasonal trends were dependent on river discharge conditions. Generally, TA was highest in the fall and winter and decreased to a minimum during spring and early summer. DIC was high and well-mixed during winter but had a strong vertical gradient in the summer due to surface productivity and subsurface respiration.

The sampling frequency employed here captured the contrast between two drought years with low river discharge (2020 and 2022) and one year with anomalously high precipitation and river discharge (2021). While late summer hypoxia dynamics in WLIS were relatively similar among all three sampling years, coincident acidification was greatest under low-flow conditions. During high-flow conditions, the DIC/TA ratio of the Housatonic River was lower, estuarine residence time decreased, and WLIS was better buffered against pH and Ω_ar decreases related to microbial respiration.

Late summer and autumn river discharge in the northeast US is expected to decrease in the coming decades, which may decrease the summertime buffering capacity of WLIS and increase estuarine residence time, leaving the estuary more vulnerable to acidification. It is unclear why a higher river discharge rate decreased riverine DIC but not TA. This motivates more detailed studies of the carbonate chemistry of the open LIS and its freshwater endmembers in relation to hypoxia dynamics as the physical and biological dynamics of the estuary change with time.