Quantifying the Effects of Nutrient Enrichment and Freshwater Mixing on Coastal Ocean Acidification

The U.S. Northeast is vulnerable to ocean and coastal acidification because of low alkalinity freshwater discharge that naturally acidifies the region, and high anthropogenic nutrient loads that lead to eutrophication in many estuaries. This study describes a combined nutrient and carbonate chemistry monitoring program in five embayments of Buzzards Bay, Massachusetts to quantify the effects of nutrient loading and freshwater discharge on aragonite saturation state (Ω). Monitoring occurred monthly from June 2015 to September 2017 with higher frequency at two embayments (Quissett and West Falmouth Harbors) and across nitrogen loading and freshwater discharge gradients. The more eutrophic stations experienced seasonal aragonite undersaturation, and at one site, nearly every measurement collected was undersaturated. We present an analytical framework to decompose variability in aragonite Ω into components driven by temperature, salinity, freshwater endmember mixing, and biogeochemical processes. We observed strong correlations between apparent oxygen utilization and the portion of aragonite Ω variation that we attribute to biogeochemistry. The regression slopes were consistent with Redfield ratios of dissolved inorganic carbon and total alkalinity to dissolved oxygen. Total nitrogen and the contribution of biogeochemical processes to aragonite Ω were highly correlated, and this relationship was used to estimate the likely effects of nitrogen loading improvements on aragonite Ω. Under nitrogen loading reduction scenarios, aragonite Ω in the most eutrophic estuaries could be raised by nearly 0.6 units, potentially increasing several stations above the critical threshold of 1. This analysis provides a quantitative framework for incorporating ocean and coastal acidification impacts into regulatory and management discussions.


Introduction
Over the past decade, ocean acidification has emerged as a major concern and negative consequence of rising atmospheric carbon dioxide (Doney, Balch, et al., 2009). Ocean acidification has been shown to significantly affect marine calcifiers, both experimentally and under natural field conditions (Barton et al., 2012;Kroeker et al., 2013Kroeker et al., , 2014Waldbusser & Salisbury, 2014). Marine calcifiers are particularly sensitive to ocean acidification's effects on calcium carbonate (CaCO 3 ) mineral saturation state (Ω), which is a measure of the thermodynamic stability of CaCO 3 . Increased carbon dioxide invasion in the surface ocean is projected to cause significant declines in global ocean pH and Ω by the end of the century , Doney, Balch, et al., 2009, Bopp et al., 2013.
Open-ocean acidification, driven by atmospheric carbon dioxide trends, is relatively well constrained for surface waters under future climate scenarios. However, projecting changes in nearshore environments is much more complicated because coastal ecosystems are affected by a number of additional processes such as freshwater discharge, eutrophication, and upwelling (Alin et al., 2012;Cai et al., 2011;Wallace et al., 2014). Given that many coastal economies rely heavily on calcifying species sensitive to ocean acidification, the future impact of ocean and coastal acidification is likely to affect coastal communities disproportionately. For instance, significant social vulnerability and economic losses are projected in the future for coastal communities that rely on shellfish as a food source or for economic revenue (Cooley et al., 2012;Cooley et al., 2015;Ekstrom et al., 2015;Mathis et al., 2015). The U.S. Northeast has been identified as particularly vulnerable both socially due to regional dependence on sensitive fisheries Hare et al., 2016) and environmentally because nearshore seawater buffering capacity is already limited due

Quantifying the Effects o Enrichment and reshwater ixing on oastal cean cidification F M C O A
Nutrien f t to low-alkalinity source water (Gledhill et al., 2015;Wang et al., 2013;Wanninkhof et al., 2015). This region is also severely impacted by eutrophication due to increasing coastal development (Valiela et al., 1997(Valiela et al., , 2016Williamson et al., 2017) and thus is vulnerable to periodic or chronic low pH and Ω Wallace et al., 2014). Some estuaries and nearshore ecosystems already experience seasonal aragonite saturation states that approach or are below 1 (Wallace et al., 2014, Gledhill et al., 2015, Wang et al., 2017, indicating environments that may be particularly stressful to calcifiers. The eutrophication problem in the U.S. Northeast is widespread and is both ecologically and economically damaging (Liu et al., 2017;Merrill et al., 2018;Nicholls & Crompton, 2018). Many coastal embayments in the U.S. Northeast have been designated impaired as pertaining to section 303(d) of the Clean Water Act (Masschusetts Department of Environmental Protection (MassDEP) 2017; Cape Cod Commission, 2017). Impaired coastal waters require the development of total maximum daily load (TMDL) limits to address and mitigate the negative effects of nutrient pollution and improve water quality. TMDL development requires observational studies combined with hydrodynamic and water quality modeling efforts to quantify target nutrient concentrations throughout the water body based on current and possible future nitrogen loading scenarios (e.g., Howes et al., 2006Howes et al., , 2012Howes et al., , 2013Howes et al., , 2014Howes et al., , 2015. To reach these TMDL limits and resulting improved nitrogen concentrations, communities develop water quality management plans that specifically target and reduce eutrophication through nitrogen loading reduction strategies (e.g., Cape Cod 208 plan, Cape Cod Commission, 2015, 2017). These management plans include provisions related to the impacts of nitrogen or phosphorous pollution but do not currently consider improvements to coastal acidification as a potential benefit to nutrient remediation. From a management perspective, addressing eutrophicationdriven coastal acidification through a regulatory framework using the Clean Water Act could benefit from direct quantification of the degree of impact of nutrient pollution on Ω (Sakashita, 2016).
Here, we describe a combined water quality and carbonate chemistry monitoring program designed to characterize the seasonal and multi-year variations in carbonate and nutrient chemistry across nitrogen loading and freshwater discharge gradients in Buzzards Bay, Massachusetts (MA). We provide an analytical framework that can be used to quantify the impacts of eutrophication on Ω and other inorganic carbon system metrics such as pH, and using the observed relationships between monitored eutrophication indicators and Ω, we project possible improvements and increases in Ω as a consequence of reductions in nitrogen loading under future loading scenarios.

Site Description
Buzzards Bay, MA is a shallow, elongated basin of approximately 11 m mean depth that is the location of an intensive, long-term citizen-science water quality monitoring program run by the Buzzards Bay Coalition, a nonprofit organization whose mission is to protect and restore the health of Buzzards Bay. This monitoring program has documented strong declines in Buzzards Bay's water quality over the past several decades that may be related to nutrient loads or to regional warming (Rheuban et al., 2016;Williamson et al., 2017).
The coastline of Buzzards Bay includes many small embayments with distinct watersheds that all receive tidal inflow from Buzzards Bay, yet have dramatically different water quality due to differences in land use and development, geomorphology, hydrodynamics, and watershed-to-estuary ratios. Land use within each embayment's watershed ranges from urban and heavily developed (e.g., New Bedford Harbor) to agriculturally dominated (e.g., Westport River) to mixed forests with sparse housing developments (e.g., Quissett Harbor; Williamson et al., 2017). In this work, we collaborated with the Buzzards Bay Coalition to augment their existing nutrient monitoring program by adding carbonate chemistry measurements.

Sample Collection and Processing
Of the approximately 200 stations regularly monitored by the Buzzards Bay Coalition, 10 stations in five embayments and a station in the open waters of Buzzards Bay were chosen for additional carbonate chemistry monitoring (Figure 1). We selected stations with the longest nutrient monitoring history, which represented endmembers in freshwater inputs and underlying nitrogen loads and water quality. All five embayments sampled also have technical reports that determine their nitrogen loads through the Massachusetts Estuaries Project (MEP, Howes et al., 2006Howes et al., , 2012Howes et al., , 2013Howes et al., , 2014Howes et al., , 2015, which uses data and sampling stations from the Buzzards Bay Coalition as validation for their nitrogen load estimates and modeling efforts.
Each embayment contained two sampling stations, denoted as "inner" or "outer," that reflected within-embayment nutrient and/or salinity gradients due to freshwater inputs, residence time, or enhanced localized nitrogen inputs. Freshwater inputs in the east branch of the Westport River (hereafter, Westport), New Bedford Harbor (hereafter, New Bedford), and the Wareham River Estuary (hereafter, Wareham) were dominated by surface waters, while West Falmouth Harbor (hereafter West Falmouth) and Quissett Harbor (hereafter, Quissett) were dominated by groundwater inputs. Estimates of nitrogen loading rates, normalized by estuary volume, varied across embayments (Table 1). Freshwater samples were also collected during July 2017 in the nontidal portion of the upper Agawam River that feeds into the Wareham River estuary complex.   Williams and Neill (2014) and Rheuban et al. (2016). Nutrient samples were processed either at the Marine Biological Laboratory or the Woods Hole Oceanographic Institution Nutrient Facility. Ancillary site characteristics (temperature, salinity, and dissolved oxygen) were measured with YSI85, YSIProPlus, or YSI600XL handheld meters calibrated prior to sampling. The full carbonate system was calculated using CO2SYS for Matlab (Lewis & Wallace 1998) using the K1 and K2 dissociation coefficients from Millero et al. (2010) and KSO 4 dissociation coefficients from Dickson (1990) with total borate from Uppstrom (1974). Apparent oxygen utilization (AOU) was calculated from dissolved oxygen measurements as the difference between saturation oxygen concentration at in situ temperature and salinity and observed oxygen concentration. Total nitrogen (TN) concentration under different loading scenarios was obtained from the MEP reports for each embayment (Howes et al., 2006(Howes et al., , 2012(Howes et al., , 2013(Howes et al., , 2014(Howes et al., , 2015.

Data Analysis
Calcium carbonate mineral saturation state, Ω, is defined using the product of calcium and carbonate ion concentrations: where K sp is the effective thermodynamic solubility product. In this analysis, we focus on Ω with respect to the carbonate mineral aragonite, since all bivalves produce an initial shell made of aragonite, and many also produce aragonite shells as adults. Hereafter, when referring to Ω, we refer to aragonite saturation state. We investigated the factors driving regional and seasonal variations in Ω computed from temperature (T), salinity (S), and carbon system variables: where the function f refers to the CO2SYS thermodynamic equations. Anomalies ΔΩ were computed for each sample relative to a reference value, Ω 0 : We used the high-salinity coastal end-member (the CBB1 site) as our reference value and determined station specific deviations from the open waters of Buzzards Bay. The mean (n = 18) water properties for the CBB1 sampling station were used for the reference salinity, temperature, DIC, and ALK and were defined as The reference saturation state is then computed from: Note that the reference value Ω 0 (1.92) differs from the mean Ω for the CBB1 (1.91) site because of the use of the reference values and nonlinearities in the thermodynamic equations. The saturation state anomalies for all samples are then estimated by a linear Taylor series decomposition (e.g., : summing the perturbation effects for multiple factors, where ∂Ω/∂X is the partial derivative of saturation state with respect to generic variable X around the reference state, and ΔX is the deviation in X from the reference value X 0 following equation (3). Here, we consider the thermodynamic effects of temperature (ΔΩ T ) and salinity (ΔΩ S ), conservative mixing of DIC and ALK of the freshwater end-member (ΔΩ mix ), and biogeochemistry factors (ΔΩ BGC ), with the addition of an error term, ε, to account for unidentified processes and nonlinearities in the carbonate system: The individual terms are approximated numerically using the CO2SYS thermodynamic code, rather than explicitly determining ∂Ω/∂X. The first two terms for the temperature and salinity perturbations arise straightforwardly from the Taylor decomposition in equation (5): The temperature and salinity terms reflect changes in Ω arising from variations in observed temperature, T, and observed salinity, S, on the thermodynamic coefficients in CO2SYS involved in K sp and calculation of [CO 3 2− ] from DIC and ALK. The salinity term also accounts for the linear variation of [Ca 2+ ] with salinity (equation (1)).
To account for correlation between ALK and DIC and strong dilution with freshwater (e.g., Fassbender et al., 2017) at inner harbor sampling stations, we use freshwater mixing (ΔΩ mix ) and biogeochemical (ΔΩ BGC ) terms rather than ΔΩ DIC and ΔΩ ALK , as has been done in a number of other studies (e.g., Doney, Lima et al., 2009 for pCO 2 ; Wang et al., 2017;Xu et al., 2017 for Ω) or a combined biophysical term (e.g., Fassbender et al., 2018 for pCO 2 ). The mixing term reflects the effects on Ω of changes of DIC and ALK as a function of salinity for two-endmember conservative mixing determined for Buzzards Bay: where DIC S and ALK S are the estimated values at the observed sample salinity S based on freshwater and CBB1 endmembers. Mixing will be linear for DIC and ALK because they are conservative tracers, but the salinity-dilution curves for derived quantities such as Ω and pH will be nonlinear. The freshwater endmember was derived from measurements collected in the upper Agawam River, which flows into the Wareham River Estuary (Figure 1) and is one of the seven main river basins flowing into Buzzards Bay. Ideally, this analysis would include measurements of DIC and ALK from multiple freshwater sources to generate mixing curves for each embayment, but the additional freshwater sampling was beyond the scope of our study. We used the mean of samples from CBB1 as the coastal seawater endmember and for measurements where salinity was beyond the mean of CBB1, we extrapolate the DIC-S and ALK-S mixing lines, but note that the extrapolation distance is relatively small. To explore the influence of the choice of endmembers on our results, we used bootstrapping techniques to resample our fresh and seawater endmembers (see Text S1 and Figures S1-S3 in the supporting information).
The biogeochemical term can be estimated as the difference between the calculated Ω at the observed sample DIC, ALK, and salinity and the Ω from the dilution curve at the sample salinity: 10.1029/2019JC015556

Journal of Geophysical Research: Oceans
where DIC, ALK, and S are the observed values. The direct estimate of four terms of equation (6) allows an explicit calculation of a residual term that arises from dropping the higher order terms in the linear Taylor decomposition (equation (5)) of a nonlinear thermodynamic system:

Overall Patterns in Carbonate Chemistry
In most embayments, we observed measurable differences in both carbonate and nutrient chemistry between inner and outer harbor stations, where the generally fresher inner harbor stations tended to have lower DIC, ALK, pH, and Ω and higher pCO 2 (Figures 2, S6-S9, and S12) and nutrient concentrations (Figures 3, S14- S17, and S19) than outer harbor stations. Carbonate system measurements showed seasonal variations due to changes in DIC, ALK, and freshwater inputs as well as temperature dependence of calculated carbonate system parameters (e.g., pH, Ω; Figure S2). Seasonal variability in DIC and ALK was evident at most sampling stations and was broadly consistent from year to year, although more years of monitoring would be necessary to define an accurate baseline climatological seasonal cycle. Buzzards Bay source water exhibited higher DIC and ALK during winter than summer, and outer harbor stations generally followed a similar pattern (Figures 2, S6, and S7). Inner harbors exhibited less consistent seasonality in DIC and ALK but nearly always had lower pH, Ω, and higher pCO 2 than outer harbor stations during summer when biological activity was high (Figures 2, S8, S9, and S12).
We observed a strong correlation between DIC and ALK across all sampling stations and seasons. The firstorder variability in these parameters was clearly related to freshwater inputs, as indicated by correlations of DIC and ALK with salinity (Figures 4 and 5). Freshwater samples from the upper Agawam River and the relatively high-salinity samples from the CBB1 station provided the endmembers for conservative mixing The observed ratio of DIC to ALK was correlated to AOU (r = −0.43, p < 0.0001) , reflecting CO 2 uptake or release and respiratory history (Figure 4). For samples with large positive AOU (i.e., O 2 uptake and CO 2 release due to both water column and benthic respiration), the DIC to ALK ratio approached 1 or was greater than 1, and thus sample Ω, [CO 3 −2 ], and pH tended to fall well below their expected values based on the mixing curve and seasonal variations in temperature (Figures 4b and 6). Samples with large negative AOU (i.e., O 2 production and CO 2 drawdown due to primary production), tended to fall above the expected mixing curves while samples with AOU close to zero tended to fall near the mixing curve ( Figure 6).
Our observations show that freshwater input is a strong driver of local carbonate chemistry, but that eutrophication also has a large influence; similar to results that have been reported by previous studies (e.g., Wallace et al., 2014, Pimenta & Grear, 2018. At many of our inner harbor stations, we observed Ω values that seasonally approached 1.0 on a regular basis, and in West Falmouth Harbor, Ω was consistently at or below 1.0 all year (Figure 2). The overall pattern was similar for calcite saturation state but shifted to higher values by approximately 1 unit (not shown). Even with the lower solubility of calcite, stations with lower salinity also experienced calcite saturation state near or less than 1 occasionally (not shown). In addition to very low Ω, we also observed lower dissolved oxygen at inner harbor stations with concentrations as low as 105.5 μmol kg −1 and 49.6 % saturation, suggesting these environments experience both acidification and dissolved oxygen stresses.

Journal of Geophysical Research: Oceans
Our sample collection was designed to target periods of high stress by sampling during the ebb tide and in the morning, potentially capturing higher freshwater discharge and less time for local production to take up DIC produced throughout the night. Higher temporal frequency measurements of pH at nearby locations suggest that local production during the day can increase pH by more than 1 unit and thus may raise Ω considerably compared to the morning, low tide conditions in this data set (see Howarth et al., 2014 for a time series from West Falmouth Harbor). Diurnal variations in pH driven by CO 2 production or consumption are usually coupled with concurrent variations in O 2 concentration that can lead not only to periods of multiple stresses (i.e., low pH, low Ω, and low O 2 at night) but also periods of beneficial conditions (i.e., high pH, high Ω, and high O 2 during the day) that may provide temporary refugia for sensitive species (Hendriks et al., 2014). However, recent experimental work on bay scallops and hard clams suggests diurnal variability in pH and O 2 simulating local production and respiration cycles does not in fact provide a refuge from acidification, and in some cases, organisms fared worse under fluctuating than chronic acidification conditions (Gobler et al., 2017).

Decomposition of Aragonite Saturation State
Developing an accurate understanding of the underlying drivers of coastal acidification is critical to determining effective management solutions. We chose to decompose Ω into four components to determine the relative contribution of temperature, salinity, freshwater mixing, and biogeochemistry on observed Ω. Figure 6b shows how these factors would affect the Ω of the Buzzards Bay study area where the reference station (CBB1) is shown by the black circle. Temperature-induced variation (ΔΩ T ) would shift Ω vertically, with lower Ω at lower T. A decrease in salinity (ΔΩ S ) from the S 0 value of CBB1 would lead to higher Ω, reflecting the offsetting effects of decreasing calcium ion concentration and K sp and increases in [CO 3 −2 ] at lower salinity due to changes in the dissociation coefficients K1 and K2 and in the activity coefficients of ions in solution. Freshwater mixing (ΔΩ mix ) would reduce Ω along the mixing curve because the freshwater endmember observed in Buzzards Bay was relatively acidic, with a DIC:ALK ratio greater than 1. Note that Ω vs. salinity mixing curves will be different in environments with high-alkalinity freshwater endmembers, such as those receiving waters from glacial runoff (e.g., Tank et al., 2016) or regions with different bedrock geology (e.g. Cai & Wang, 1998). Local biogeochemistry (ΔΩ BGC ) would shift Ω vertically with no change in salinity, with the sign and magnitude of this change dependent on the stoichiometry of the metabolic processes (e.g., Soetaert et al., 2007). Each sample's overall perturbation from the coastal source water (ΔΩ Total ) is assumed to be a combination of these four factors and a residual term (ε) that includes any unaccounted for processes and limitations of linear approximations of nonlinear systems.
The term ΔΩ Total varied both between embayments and within each embayment, and inner harbor stations typically had larger negative ΔΩ Total than outer harbor stations (Figure 7). Temperature and salinity had a relatively small overall impact on the observed range of ΔΩ Total , a result consistent with similar analyses of Mid-Atlantic Bight shelf/slope waters (Xu et al., 2017). However, at outer harbor stations where Buzzards Bay water had a large influence on the local carbonate system but the overall variation in Ω was small, seasonal variability in temperature and salinity caused ΔΩ T and ΔΩ S to, in a few cases, be of a similar magnitude as ΔΩ Total (e.g., outer Quissett). Freshwater mixing had a clear and very strong control on the Figure 6. Top panels: aragonite saturation state (Ω) vs. salinity, and conceptual diagram illustrating the decomposition of Ω as a function of salinity. The solid line for all panels is calculated from the conservative mixing line for DIC and ALK using the dataset mean temperature, while dashed lines are calculated using the maximum and minimum observed temperatures. The dot indicates the CBB1 reference station. Each perturbation process is shown as ΔΩ x , where ΔΩ x describes the perturbation of Ω due to temperature (T), salinity, (S), conservative freshwater mixing (mix), and biogeochemical processes (BGC) (equations (6)-(10)). Arrows indicate the expected direction caused by each perturbation process. Bottom panels: carbonate ion concentration (CO 3 −2 ), and pH vs. salinity with the lines calculated as above from the conservative mixing lines. All points are colored by apparent oxygen utilization (AOU). Data points are given unique symbols that represent specific embayments that are consistent in the figures throughout the manuscript.

Journal of Geophysical Research: Oceans
variations in Ω where ΔΩ mix was consistent with differences in salinity and expected freshwater discharge across embayments. Stations with high surface or groundwater discharge tended to have relatively low salinity (Westport, Wareham, and West Falmouth inner harbor stations), and thus larger ΔΩ mix (Figure 7). The residual term was correlated to salinity such that the magnitude was larger for fresher measurements, but for the majority of our samples, the residual term was small relative to the magnitude of ΔΩ Total (Figures 7 and S5).
In most embayments, ΔΩ BGC was a large component of ΔΩ Total . Outer harbor stations tended to have smaller ΔΩ BGC than inner harbor stations, due to carbonate system chemistry that was similar to the open waters of Buzzards Bay, yielding smaller deviations from the estimated mixing curves. ΔΩ BGC is likely driven by a combination of the relative magnitudes of metabolic fluxes and the local water residence times. Small metabolic fluxes with short residence times would yield small ΔΩ BGC , while increasing metabolic fluxes in the same location would increase ΔΩ BGC . Small flux rates in locations with long residence times may also lead to increased ΔΩ BGC , and we would expect large flux rates in areas of long residence time to lead to the largest ΔΩ BGC. We observed correlations between station AOU and ΔΩ BGC (r = −0.46, p < 0.0001, all data; r = −0.57, p < 0.0001, summer, Figure 8), consistent with biological processes being the driver of the ΔΩ BGC component, especially during periods of warmer water temperatures and likely higher biological activity. The slope of the observations during warm periods (temperature > 15°C, −0.0056 ± 0.0007 ΔΩ BGC /AOU, (1/(μmol kg −1 )) was also consistent with aerobic respiration and production using a Redfield ratio of DIC and ALK to O 2 (−0.0097 ΔΩ BGC /AOU, (1/(μmol kg −1 ); Figure 8).
Although summertime ΔΩ BGC was consistent with oxygen production and aerobic respiration there was still considerable variation, especially in samples collected during cooler temperatures, that could be driven by other processes. First, our dilution curves do not capture embayment-specific or seasonal variability in endmember composition. In particular, it is possible that the carbonate chemistry of freshwater endmembers for  (6)) as: the total perturbation from the CBB1 reference station (Total), the perturbations due to temperature (T), salinity (S), freshwater mixing (Mix), and biogeochemical processes (BGC), and the estimated residual term (Residual). Bottom panels (blue) are from outer harbors, and top panels (red) are from inner harbor stations, the far left panels in both the top and bottom rows are from the CBB1 reference station. Dashed lines indicate the median ΔΩ Total for each sampling station.

10.1029/2019JC015556
Journal of Geophysical Research: Oceans embayments across Buzzards Bay differs substantially from the Agawam River source water used here or that there may also be seasonal changes in surface or groundwater endmember DIC or ALK values. Both of these situations could lead to changes in the relative magnitudes of ΔΩ mix or ΔΩ BGC by altering the freshwater endmember of the dilution curve. The freshwater endmember DIC and ALK values from the upper Agawam River (517 ± 131 and 409 ± 103 for DIC and ALK, respectively, n = 6) agreed well with riverine carbonate chemistry data from the GLORICH database (Hartmann et al., 2019, see Text S3) averaged from a 100 km radius surrounding Buzzards Bay (665 ± 279 and 473 ± 225 [SD, n = 16 stations] for DIC and ALK, respectively). A bootstrapping error analysis found that the ΔΩ BGC term was not sensitive to varying the mean freshwater and seawater endmember DIC and ALK within the observed range (e.g., ± 0.045 Ω units, see Text S1, Figures S1-S3). Additionally, for simplicity in this analysis, we used a single seawater endmember calculated from the mean DIC, ALK, S, and T measured in the open waters of Buzzards Bay. We did observe seasonality in the Buzzards Bay (CBB1) data and outer harbor measurements (Figures 2, S6-S9, and S12), which could indicate seasonal variation in either the coastal water masses or biogeochemical fluxes that inherently alter the carbonate chemistry of the Buzzards Bay source seawater. However, these observed seasonal variations are relatively small, and reanalyzing the dataset with a seasonally varying seawater endmember did not substantially change the results of this analysis (data not shown). Furthermore, we do not believe we have enough data from individual embayments to generate seasonal or embayment-specific dilution curves.
Second, anaerobic biological processes such as denitrification or sulfate reduction, nonbiological redox reactions, or significant calcification or dissolution may alter the seawater DIC:ALK ratio Glud, 2008;Liu et al., 2017;Soetaert et al., 2007) and thus lead to the variability observed in the ΔΩ BGC to AOU relationship. These processes could occur in the water column or sediments, which often represent a large portion of the total ecosystem fluxes in shallow coastal environments (e.g., Long et al., 2019). The bulk water column measurements collected for this study integrate the effects of both water column and sedimentary processes but do not allow us to estimate the relative contribution of sediments and water column processes to the ΔΩ BGC term. Although we highlight the correlation between ΔΩ BGC and AOU as an indication of one possible biological process affecting the ΔΩ BGC term, the analysis described below does not rely on the relationship between ΔΩ BGC and AOU relationship in any way nor do we make any assumptions on the underlying process driving the observed patterns.

Influence of Eutrophication on Aragonite Saturation State
Regardless of the mechanism, we hypothesize that of the set of factors we use to describe Ω, the ΔΩ BGC term is most likely influenced by eutrophication. We observed a strong correlation between July and August station mean ΔΩ BGC and mean TN concentration (r = −0.82, p = 0.002, Figure 9). Since TN concentration reflects eutrophication status in coastal Buzzards Bay (Rheuban et al., 2016), this is consistent with a strong eutrophication influence on ΔΩ BGC . Annual mean ΔΩ BGC was also strongly correlated to annual mean TN (r = −0.85, p = 0.001, not shown), but we focus our analysis on summer mean ΔΩ BGC and TN because historical measurements of water quality in Buzzards Bay were typically only collected during the summer, and regional water quality  The pattern in both TN and ΔΩ BGC followed expected nitrogen loading gradients ( Figure 9); bays with high nitrogen loading and resultant high rates of respiration and organic matter decomposition had strongly negative ΔΩ BGC values. Within-embayment variation (inner vs. outer harbor ΔΩ BGC ) was also consistent with this mechanism, with inner harbors tending to fall further from the Buzzards Bay reference station ( Figure 9). The inner station in New Bedford was the exception to this relationship, in part because we observed very different chemistry in summer 2015 and summer 2016. In 2016, DIC and ALK at the inner station were much lower than expected based on the observed salinity. The 2016 DIC:ALK ratio was also lower than in 2015, leading to Ω that was considerably higher (2.54 ± 0.29, summer 2016) than observed during other sampling periods in New Bedford (1.35 ± 0.28 for summer 2015, 1.57 ± 0.42 across all data). The inner New Bedford mean from 2015 agreed well with other measurements but without an understanding of the cause of the summer 2016 persistent high Ω, we hesitate to exclude it from our analysis. This high degree of variability between the two summers sampled in New Bedford illustrates the importance of long-term monitoring of the coastal carbonate and nutrient systems in order to fully develop baseline climatologies that will allow for the exploration of potential future change.
Regional water quality managers have very little control over temperature, salinity, and freshwater mixing; however, state and local managers can use the regulatory framework of the Clean Water Act to reduce nitrogen pollution and the negative consequences of eutrophication. As part of statewide efforts to improve coastal water quality, the MEP (https://www. mass.gov/guides/the-massachusetts-estuaries-project-and-reports) was designed to evaluate estuarine health and develop TMDLs for approximately 70 embayments across the Massachusetts coastline, including the five embayments sampled as part of this study (Howes et al., 2006(Howes et al., , 2012(Howes et al., , 2013(Howes et al., , 2014(Howes et al., , 2015. These MEP reports modeled TN concentration under proposed future nitrogen loading scenarios at most of our sampling stations because the MEP process used the same Buzzards Bay Coalition historical data and sampling stations as validation for their water quality modeling efforts. Using modeled TN concentrations, the relationship observed between ΔΩ BGC and TN (Figure 9) can then provide an estimate of the impacts of future nitrogen loading reduction scenarios on station Ω. Of the five harbors sampled, four of the inner stations show projected increases in Ω under both the reduced nitrogen loading TMDL and the "no anthropogenic Figure 10. Average aragonite saturation state (Ω) from observations (blue) and estimated based on two nitrogen loading reduction scenarios: the total maximum daily load (TMDL, red) and no anthropogenic nitrogen loading (yellow). TN concentration was obtained from reports for each harbor generated through the Massachusetts Estuaries Project. TN was determined through water quality modeling based on nitrogen loads at proposed TMDLs and no anthropogenic nitrogen loading scenarios. Dashed line illustrates Ω = 1. Error bars are standard error, estimated by propagating the variability in observed Ω, the uncertainty associated with the DIC and ALK mixing curves, and the uncertainty in the relationship between total nitrogen concentration (TN) and ΔΩ BGC .  Howes et al. (2006Howes et al. ( , 2012Howes et al. ( , 2013Howes et al. ( , 2014Howes et al. ( , 2015. Error shown is standard error (n = 12, 4, 4, 4, 12, respectively). ΔΩ BGC for the TMDL and No Load scenarios is determined from the TN concentration and the relationship in Figure 10. Ω is then estimated from equations (6) and (3).

10.1029/2019JC015556
Journal of Geophysical Research: Oceans load" scenarios, with the exception of New Bedford inner harbor, which shows declines that are not statistically different than the observations ( Figure 10, Table 2). Harbors with better water quality, higher observed Ω, and lower present-day estuary-volume-normalized nutrient loads see smaller or little improvement. We estimate that Ω in morning, low-tide conditions in the most eutrophic estuaries may increase (i.e., improve) by nearly 0.5 units under a TMDL scenario, and more than 0.6 units under a no anthropogenic load scenario ( Figure 10, Table 2).
These projected increases in station summer mean Ω could lead to dramatic improvements in West Falmouth and Wareham inner harbors in particular, where Ω may increase from conditions close to or well below an Ω of 1 to above the critical thermodynamic equilibrium (Figure 10). Our sampling program was designed to target the periods with the greatest stresses (morning, low-tide), and as it is possible that our seawater endmember is also influenced by eutrophication, these improvements would represent a conservative lower bound for average Ω in these bays. Model estimates for the northeast shelf suggest that Ω has declined by more than 0.5 units since the industrial revolution , Environmental Protection Agency (EPA), 2016. A simple calculation of changes to Ω for Buzzards Bay from preindustrial conditions (see Text S2) illustrates that the Ω improvements (e.g., increases) that could be achieved through reduced nitrogen loading at the most eutrophic stations are likely to be larger than the Ω declines associated with past or near-future atmospheric CO 2 increases (Table S1). Many shellfish species important to the culture and economy of the U.S. Northeast, e.g., hard clam, eastern oyster, and bay scallop, have shown sensitivity to Ω (see summary in Gledhill et al., 2015) and most laboratory studies show that lowered Ω values lead to reduced growth, calcification, or survival rates even for experimental conditions with Ω larger than 1. Thus, these potential improvements to Ω could facilitate improved shellfish survival and growth even for embayments not yet undersaturated with respect to aragonite.

Conclusions
Application of management strategies to improve coastal acidification is complex, as there are many drivers that determine local carbonate chemistry. Furthermore, developing regulatory loading limits through the Clean Water Act for dissolved carbon in coastal environments is not likely a suitable strategy to improve coastal acidification. However, the connection between eutrophication and coastal acidification is clear (Cai et al., 2011;Howarth et al., 2011;Wallace et al., 2014;Pimenta & Grear 2018). As a result, more immediate and larger improvements to the coastal acidification problem can likely be gained by reducing eutrophication and the subsequent organic matter respiration and CO 2 release. Additionally, in highly eutrophic estuaries, the increases to Ω gained though nutrient management could be larger than the estimated decline associated with a change from preindustrial atmospheric CO 2 to present day conditions and up to nine fold larger than those associated with Ω declines expected over the next decade due to atmospheric CO 2 trends (Table S1). The quantitative framework described in this study allowed us to estimate the relationship between proxies for nitrogen loading and coastal acidification in order to predict the possible improvements in coastal carbonate chemistry following nitrogen loading reduction scenarios. We envision this framework will be useful for local and regional decision making to allow quantitative metrics of potential improvements in coastal acidification to be incorporated in the discussion of nutrient impairment and TMDL attainment strategies.