Decadal-Scale Acidification Trends in Adjacent North Carolina Estuaries: Competing Role of Anthropogenic CO2 and Riverine Alkalinity Loads

Decadal-scale pH trends for the open ocean are largely monotonic and controlled by anthropogenic CO₂ invasion. In estuaries, though, such long-term pH trends are often obscured by a variety of other factors, including changes in net metabolism, temperature, estuarine mixing, and riverine hydrogeochemistry. In this study, we mine an extensive biogeochemical database in two North Carolina estuaries, the Neuse River estuary (NeuseRE) and New River estuary (NewRE), in an effort to deconvolute decadal-scale trends in pH and associated processes. By applying a Generalized Additive Mixed Model (GAMM), we show that temporal changes in NewRE pH were insignificant, while pH decreased significantly throughout much of the NeuseRE. In both estuaries, variations in pH were accompanied by increasing river discharge, and were independent of rising temperature. Decreases in bottom-water pH in the NeuseRE coincided with elevated primary production in surface waters, highlighting the importance of eutrophication on long-term acidification trends. Next, we used a simple mixing model to illustrate the impact of changing river discharge on estuarine carbonate chemistry. We found that increased riverine alkalinity loads to the NewRE likely buffered the impact of CO₂-intrusion-induced acidification. In the NeuseRE, however, elevated dissolved inorganic carbon loads further decreased the buffering capacity, exacerbating the effects of CO₂-intrusion-driven acidification. Taken together, the findings of this study show that future trajectories in estuarine pH will be shaped by complex interactions among global-scale changes in climate, regional-scale changes in precipitation patterns, and local-scale changes in estuarine biogeochemistry.

Introduction

Anthropogenic activities have caused dramatic increases in atmospheric CO₂ levels, from a pre-industrial level of <300 ppm to the current value of 407 ppm (Keeling and Keeling Charles, 2017). Much of the open ocean is under-saturated in CO₂ with respect to the atmosphere, causing a net uptake of CO₂, mitigating approximately 30% of anthropogenic emissions (Le Quéré et al., 2016). This absorption of atmospheric CO₂ is leading to the increase in aqueous CO₂ concentration {[CO_2(aq)]}, which affects the speciation of inorganic carbon in the ocean, shifting the following reaction (Eq. 1) toward the right, favoring the formation of H⁺ and lowering the pH. Accordingly, anthropogenic CO₂ emissions have caused ocean pH to fall at a rate of approximately 0.001–0.0025 yr^-1 (Doney et al., 2009; Bates et al., 2014; Takahashi et al., 2014).

Excess H⁺ produced by the dissolution of CO_2(aq) pairs with CO₃^2- to form HCO₃^-, resulting in a net decrease in [CO₃^2-] and saturation states (Ω) of various carbonate minerals, adding stress to marine calcifying organisms. This process is collectively known as ocean acidification, and may directly damage calcifying organisms, further decrease the fitness of commercially valuable groups, by directly damaging shells or compromising early development and survival.

Over very long time scales (i.e., >100,000 years), it is expected that increased continental weathering and subsequent delivery of alkalinity via rivers will partially mitigate the coastal and open ocean acidification signal (Doney and Schimel, 2007; Müller et al., 2016). Over shorter time scales, though, riverine inputs of freshwater and related solutes may cause varied effects on the carbonate system in estuaries, resulting in contrasting pH trajectories (Aufdenkampe et al., 2011; Hu and Cai, 2013; Hu et al., 2015; Müller et al., 2016). For example, changing local water budgets resulted in a range of temporal trajectories in riverine alkalinity loading for a set of Northwest Gulf of Mexico estuaries, where estuaries experiencing decreasing alkalinity loads also became acidified (Hu et al., 2015). Similarly, long-term alkalinity increases in the Baltic sea have partially or wholly counteracted anthropogenic CO₂-induced acidification in that system (Müller et al., 2016). However, the impact of short-term variations in river discharge on the estuarine carbonate system is less clear (Mote et al., 2008; Phlips et al., 2011; Couldrey et al., 2016; Baumann and Smith, 2017; Siam and Eltahir, 2017).

A variety of biogeochemical processes will interact with these watershed-scale drivers, either enhancing or counteracting acidification trends in estuaries. For example, nutrient-enhanced primary production has been linked with elevated pH in well-mixed mesocosm experiments (Nixon et al., 2015), but net production and respiration are decoupled in stratified estuaries, where the latter process decreases bottom water pH (Feely et al., 2010; Cai et al., 2011; Waldbusser and Salisbury, 2014). Still, net ecosystem metabolism in productive coastal systems like seagrass beds has been suggested to lessen long-term coastal ocean acidification (Unsworth et al., 2012), but increase the short-term extremes in pH (Cyronak et al., 2018; Pacella et al., 2018). H₂S oxidation and subsequent H₂SO₄ additions may exacerbate acidification near the pycnocline of stratified systems (Cai et al., 2017). Anaerobic production of CO₂ and alkalinity in wetlands can act to simultaneously alkalize and acidify pore-water, leading to variable impacts on estuarine pH (Hu and Cai, 2011; Wang Z.A. et al., 2016). At the same time, estuarine waters are mixing with oceanic water that is also becoming enriched in CO₂ due to anthropogenic CO₂ emissions.

The above natural and anthropogenic processes will combine non-linearly to affect pH, to an extent that may be exacerbated by the typically poor buffering capacity of estuarine water (Feely et al., 2010; Cai et al., 2011; Hu and Cai, 2013; Jeffrey et al., 2016; Van Dam et al., 2018a). Future trajectories in estuarine acidification will depend on the interactions among a variety of complex drivers, including the invasion of anthropogenic CO₂, changes in freshwater delivery, and changes in net ecosystem metabolism. However, we do not yet know how these factors will interact over decadal time scales, or across different estuaries. The estuaries of eastern North Carolina, United States, offer a unique opportunity to compare the impacts of these factors on decadal-scale acidification trends. Both New and Neuse River estuaries, NC, experience the same climate, and have seen increases in river discharge that are of a similar magnitude over the past 10 years, but contrasting watersheds properties cause dissolved inorganic carbon (DIC) and total alkalinity (TA) loading to differ (Van Dam et al., 2018a). In this study, we investigate decadal trends in pH and associated parameters in these two adjacent estuaries, then use a simple modeling approach to examine the sensitivity of estuarine acidification to above three factors.

Materials and Methods

Site Description

The Neuse River Estuary (NeuseRE) and New River Estuary (NewRE) are microtidal (tidal amplitude < 1 m), lagoonal systems located in the eastern North Carolina coastal plain (Figure 1), both of which have relatively long fresh water residence times [average 58 (±34) and 46 (±27) days, respectively (Van Dam et al., 2018a)]. Both estuaries are shallow, with average water depths of 1.8 and 2.7 m in the NewRE and NeuseRE, respectively. While the NewRE is connected directly with the coastal ocean through the New River inlet, a string of barrier islands (Outer Banks, NC) with narrow tidal inlets restricts water exchange between the NewRE and the ocean. The long residence time and limited tidal exchange has interacted with a history of anthropogenic nutrient enrichment to drive both estuaries toward eutrophy. Historically, average primary productivity (PPR) has been approximately 200–500 g C m^-2 yr^-1 in the NeuseRE (Boyer et al., 1993; Mallin et al., 1993) and <250 gC m^-2 yr^-1 in the NewRE (Mallin et al., 2005). A key difference between these two microtidal, lagoonal estuaries is their size; the surface area of the NeuseRE (352 km²) is approximately 5 times that of the NewRE (79 km²), while the Neuse River watershed (15,700 km²) is over 15 times larger than that of the New River (1,024 km²). In contrast with the heavily agricultural Neuse River watershed, much of the New River watershed contains impervious surfaces, causing river discharge in this system to be relatively ‘flashy’ in response to storm events (Hall et al., 2012; Peierls et al., 2012). Submarine groundwater discharge is a very small contributor to total freshwater loads in both NewRE (Crosswell et al., 2017) and NeuseRE (Fear et al., 2007; Null et al., 2011).

FIGURE 1

Figure 1. Site map, including sampling stations for both estuaries, and watershed extents for both rivers.

Long-Term Data Collection

Long-term data were acquired from two monitoring programs established in each estuary from 2005 to 2017. Associated with the Defense Coastal Estuarine Research Program (DCERP), monthly sample collection began in the NewRE in the Fall of 2007. The goal of this project was to assess the impacts of land use in the watershed on water quality in the NewRE. Likewise, bi-weekly to monthly sample collection occurred in the NeuseRE beginning in 1999, as a part of the Neuse River Estuary Modeling and Monitoring Project (ModMon). ModMon is a collaborative effort of both universities and the state of North Carolina, aiming at the collection of water-quality compliance data. For DCERP and ModMon programs, surface (∼0.2 m) and bottom-water (∼0.5 m above bottom) samples were collected at each station using a diaphragm pump, and were immediately stored in a cooler. For each surface sample, an estimate of phytoplankton PPR was made using the ¹⁴C uptake method, and in-vitro Chlorophyll-a fluorescence (Chl-a IV) was determined using acetone extraction (Paerl et al., 1998; Hall et al., 2012). The ¹⁴C uptake method generates PPR estimates that lie somewhere between gross and net PPR. This is because a fraction of the ¹⁴C incorporated by phytoplankton is respired during the incubation period, and lost as CO₂ before the final scintillation count. However, this method has remained stable over time in both ModMon and DCERP programs, allowing us to reliably assess trends in PPR over time.

At each DCERP and ModMon station, vertical profiles of pH, dissolved oxygen (DO), salinity (Sal), and temperature (Temp) were collected at a vertical resolution of 0.5 m (YSI 6600 multi-parameter sonde, Yellow Springs Inc, Yellow Springs, OH). The pH sensor (YSI model 6561) was calibrated against NIST-traceable NBS buffers, and has a manufacturer’s suggested accuracy of ±0.2 units, which is far less than the typical spatial variability in pH (Van Dam et al., 2018a). Once water samples were returned from the field, a small aliquot was immediately stored, unpreserved, in a 20 mL scintillation vial with no headspace, and was analyzed for DIC within 24 h using a Shimadzu TOC-5000A in inorganic carbon mode. Because DIC in these unpreserved water samples is often lower than when preserved with HgCl₂, or by filtration, we applied a correction factor (Crosswell et al., 2012) to DIC values obtained from the long-term monitoring dataset (DIC = 7.18 + 1.16 × DIC_unpreserved). We then calculated TA from this corrected DIC value, pH, Temp, and Sal using CO2SYS (Lewis and Wallace, 1998), with the carbonic acid dissociation constants of Millero (2010), and the NBS scale for pH.

Statistical Approach

We calculated trends in pH/PPR/Chl-a using a Generalized Additive Mixed Modeling (GAMM, Eq. 2), a method which has been previously applied to calculate multi-decadal fCO₂ trend in various coastal environments (Wang H. et al., 2016, 2017; Reimer et al., 2017). Briefly, a penalized spline was first used to fit the seasonal cycle, and a harmonic function was adopted to fit daily cycles, then linear regressions were used to model the variability resulting from salinity and DO change. Sampling date was also included in GAMM as a linear effect, and its coefficient represents the pH/PPR/Chl-a long-term change. The method also weighted the observations using explicit models for heteroscedasticity to account for an observed unequal variance in salinity.

$\mathrm{p}\mathrm{H}\mathrm{/}\mathrm{P}\mathrm{P}\mathrm{R}\mathrm{/}\mathrm{C}\mathrm{h}\mathrm{l}\mathrm{-}\mathrm{a}\mathrm{=}\mathit{f}\mathrm{\left(}\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{\right)}\mathrm{+}\mathit{f}\mathrm{\left(}\mathrm{d}\mathrm{a}\mathrm{i}\mathrm{l}\mathrm{y}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{\right)}\mathrm{+}\mathit{f}\mathrm{(}\mathrm{D}\mathrm{O}\mathrm{,}\mathrm{S}\mathrm{a}\mathrm{l}\mathrm{)}\mathrm{+}\mathit{f}\mathrm{(}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}\mathrm{-}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{d}\mathrm{)}\mathrm{+}\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{}\mathrm{\left(}\mathrm{2}\mathrm{\right)}$

Note, the “best” model was selected based on smallest AICc values (Supplementary Table S1 and Supplementary Figure S1). Trends were assessed for the time period 2005–2017 in the NeuseRE, and for 2007–2017 in the NewRE. Salinity and temperature trends were calculated after removing the seasonal cycle only. Pre-2005 data were not used for the NeuseRE because: (1) residuals were not normally distributed before 2005, (2) different pH probe was used before 2005, and data were more scattered, (3) surveys were conducted in early morning, causing pH measurements to be biased toward lower values (based on typical diel variations in these estuaries), and (4) data do not exist for the NewRE before 2007.

Box Model

In this study, we used a simple box model approach to address the relative importance of the following three factors on long term estuarine acidification trends: (1) continuous anthropogenic CO₂ invasion, (2) varying PPR, and (3) changing river discharge (and TA/DIC loads). This model is intended to represent average conditions over yearly to decadal-time scales, rather than shorter-term variations in pH, which can be very difficult to predict. As future changes in river discharge are very uncertain, we estimated a range of river discharge (Q_R) by estimating the average discharge (Q_R) from Supplementary Table S2 over time, extrapolating forward 70 years, and backward 20 years. This range of Q_R (New River: 1.3–8.9 m³ s^-1, Neuse River: 82 – 221 m³ s^-1) is intended to include current conditions, and what could be expected for flood and drought years in the next century. The riverine end-member concentrations of DIC (DIC_R) was set based on the empirical relationships with Q_R determined in Van Dam et al., 2018b (Supplementary Table S2). We set the riverine TA (TA_R) value based on the linear relationship between measured DIC and calculated TA for the station closest to the river mouth, when salinity was <5 for the NewRE, or <2 for the NeuseRE. Those equations are also shown in Supplementary Table S2. The r² for these relationships were 0.96 for the NeuseRE, and 0.97 for the NewRE. We assumed a constant oceanic TA end-member (TA_O) of 2500 μmol kg^-1, and increasing oceanic DIC (DIC_O) as a function of increasing atmospheric CO₂ (i.e., DIC_O ∼pCO₂, TA). These end-member DIC and TA values were mixed within the estuary, allowing us to calculate pH under different atmospheric CO₂ levels.

After setting the end-member chemistry, we assumed that the main factors affecting DIC concentration in surface and bottom water was net ecosystem metabolism. The NewRE and NeuseRE are both stratified estuaries, where net ecosystem production and respiration are vertically decoupled, meaning that some fraction of net production in the surface is exported to bottom waters where it is respired. Hence, we used average net ecosystem productivity (, mol m^-2 yr^-1) from Herrmann et al. (2015), and assumed that a certain fraction of net production was respired in bottom waters (R). Note that the sign of was reversed from that given in Herrmann et al. (2015), such that positive values indicate net heterotrophy. Given the paucity of data regarding bottom-water respiration in these systems, we set R at 40% (R = 0.4 ). Variations in R do not affect the general trend of modeled pH with variations in Q_R or pCO₂, but does cause moderate changes in the absolute value of modeled pH. As an example of these variations, modeled pH results are shown for R values of 0.25 (a), 0.5 (b), and 0.75 (c) in Supplementary Figure S4.

Areal was converted into a concentration change by the following equation where z_ave is the average water depth, τ_FW is the residence time (years), and F_mix is the mixing fraction, given by the ratio of Sal_E/Sal_O (i.e., F_mix = 10/35 = 0.29 for oligohaline region). Multiplying by F_mix scales the DIC loss due to air-water exchange from the whole estuary down to the region of interest.

${\mathrm{N}\mathrm{E}\mathrm{P}\mathrm{E}}_{\mathrm{E}}\mathrm{\left(}\mathrm{\mu }\mathrm{M}\mathrm{\right)}\mathrm{=}{\mathrm{10}}^{\mathrm{3}}\mathrm{\times }\mathrm{\left(}\frac{\overline{\mathrm{N}\mathrm{E}\mathrm{P}}\mathrm{\times }\mathrm{F}\mathrm{W}\mathrm{\times }{\mathrm{F}}_{\mathrm{M}\mathrm{i}\mathrm{x}}}{{\mathrm{z}}_{\mathrm{a}\mathrm{v}\mathrm{e}}}\mathrm{\right)}\mathrm{}\mathrm{\left(}\mathrm{3}\mathrm{\right)}$

With DIC consumption/production terms defined, surface and bottom water DIC in the estuary (DIC_S and DIC_B respectively) was calculated assuming conservative mixing, with the addition of net metabolism and air-sea exchange:

${\mathrm{D}\mathrm{I}\mathrm{C}}_{\mathrm{s}}\mathrm{=}\mathrm{(}{\mathrm{D}\mathrm{I}\mathrm{C}}_{\mathrm{O}}\mathrm{\times }{\mathrm{F}}_{\mathrm{M}\mathrm{i}\mathrm{x}}\mathrm{)}\mathrm{+}\mathrm{(}{\mathrm{D}\mathrm{I}\mathrm{C}}_{\mathrm{R}}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{F}}_{\mathrm{M}\mathrm{i}\mathrm{x}}\mathrm{)}\mathrm{)}\mathrm{-}{\mathrm{N}\mathrm{E}\mathrm{P}}_{\mathrm{E}}\mathrm{}\mathrm{\left(}\mathrm{4}\mathrm{\right)}$

$\mathrm{D}\mathrm{I}\mathrm{C}\mathrm{=}\mathrm{(}{\mathrm{D}\mathrm{I}\mathrm{C}}_{\mathrm{O}}\mathrm{\times }{\mathrm{F}}_{\mathrm{M}\mathrm{i}\mathrm{x}}\mathrm{)}\mathrm{+}\mathrm{(}{\mathrm{D}\mathrm{I}\mathrm{C}}_{\mathrm{R}}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{F}}_{\mathrm{M}\mathrm{i}\mathrm{x}}\mathrm{)}\mathrm{)}\mathrm{-}\mathrm{\left(}\mathrm{R}\mathrm{\right)}\mathrm{}\mathrm{\left(}\mathrm{5}\mathrm{\right)}$

Total alkalinity in the estuary (TA_E) was determined by conservative mixing of river and ocean end-members. We also accounted for the small amount of TA production (consumption) that occurs during the production (remineralization) of Redfield organic matter, by assuming TA changes with DIC at a ratio of -17/106.

From these modeled TA and DIC values, estuarine pH (pH_est) was calculated in CO2SYS using the carbonic acid dissociation constants of Millero (2010), at a constant temperature of 20°C, very close to the median temperature in the NeuseRE (20.5°C) and NewRE (20.8°C). A constant output temperature was used to control for the effect of changing temperature on pH and pCO₂.

Results and Discussion

Physical Condition Change

From 1996 to 2017, Neuse river discharge at Ft Barnwell increased at an average rate of 1.55 m³ s^-2 yr^-1, or 2.8% yr^-1, while New River discharge at Gum Branch increased at 0.084 m³ s^-2 yr^-1 (1.4% yr^-1) between 2007 and 2017. Variations in monthly averaged discharge during the Spring and Winter were significantly correlated with El Nino Southern Oscillation (ENSO) patterns, such that positive ENSO phase was associated with increases in discharge (Supplementary Figure S2). However, ENSO phase could not explain changes in discharge during the Fall and Summer. The observed increases in river discharge were reflected in a large and significant freshening trend in both estuaries. Averaged across all surface stations, salinity decreased at a rate of 0.002 to 0.35 yr^-1 in the NeuseRE, and 0.18 to 0.96 yr^-1 in the NewRE (Figure 2D). In both estuaries, trends in bottom water salinity followed those from surface water.

FIGURE 2

Figure 2. Post-2005 trends for pH (A), PPR (B), Chl-a (C), Salinity (D), and Temperature (E) for surface (+) and bottom water (∇). Trends are shown relative to distance downstream, which is the distance from the head of the estuary divided by the total length of the estuary. Colors indicate the p-value for a t-test on the slope of the trend line (green: p < 0.01; yellow: 0.01 < p < 0.1; red: p > 0.1).

Despite generally non-significant changes in Chl-a (Figure 2C), PPR in the middle and lower NeuseRE increased significantly from 2005 to 2017, at a rate of approximately 30 mg C m^-3 hr^-1 decade^-1. Elevated river discharge has been shown to drive a seaward shift in the location of peak phytoplankton biomass (Peierls et al., 2012), so it is not surprising that increased river inputs were accompanied with elevated PPR in the lower NeuseRE (Figure 2B). Coinciding with this increasing surface water PPR in the NeuseRE, bottom water pH fell significantly at approximately ∼0.02 yr^-1 (0.2 decade^-1), while surface water pH trends were generally non-significant. Trends in bottom water DO were not significantly different from zero in the NeuseRE (Supplementary Figure S3), indicating that the observed decreases in pH may have been accelerated by factors beyond aerobic decomposition.

Trends in both PPR and pH for the NewRE were not significant (Figures 2A,B). Meanwhile, changes in temperature over the study period were generally insignificant and variable, between -0.09 and 0.28°C yr^-1 (average = 0.03 ± 0.06°C yr^-1), similar to the trends observed in the coastal ocean (Reimer et al., 2017). Because pH decreases with temperature (holding DIC and TA constant), it is possible to assess the thermal effect on pH changes. At a TA and DIC of 900 and 1000 μmol kg^-1 respectively, an average temperature increase of 0.03°C yr^-1 results in a calculated pH decrease of -0.00024 yr^-1. This is approximately ∼2 orders of magnitude below the observed trends in pH (Figure 2A), indicating that the effect of warming waters on acidification was negligible, relative to other factors.

Similar to the observed trend of acidification in the NeuseRE, which was accompanied by consistent DO, pH in the NewRE was stable (Figure 2A), while bottom water DO in the upper NewRE showed significant decreases (Supplementary Figure S3). This discrepancy between long term trends in DO and pH may be partially explained by the difference in air-water exchange for O₂ and CO₂, which is more rapid for O₂ than for CO₂. These counterintuitive trends in long-term acidification may be explained by changes in chemical buffering that served to accelerate acidification in the NeuseRE and to buffer acidification stress in the NewRE. Further explanation for this is given in the following section.

Box Model Results

In both estuaries, a doubling of oceanic pCO₂ from 400 to 800 μatm (reflecting estimated increases in atmospheric CO₂ over the next century) caused modeled surface water pH_est to fall by 0.094–0.14 in the NewRE, and by 0.037–0.23 in the NeuseRE. This modeled future decrease in pH over the next century corresponds to a rate of change of 0.004–0.023 decade^-1, which is within the range expected for the open ocean, but is an order of magnitude slower than the observed trend for the NeuseRE (Figure 3). This suggests that future trends in acidification will be driven by multiple interacting factors including the invasion of anthropogenic CO₂, and other factors including changing riverine TA loads or NEM. In the upper NeuseRE, elevated discharge was associated with decreasing surface water pH_est (pH_est at high Q_R – pH_est at low Q_R = ΔpH_est = -0.086). Throughout the rest of the NeuseRE, increases in Q_R corresponded with increasing surface pH_est, such that ΔpH_est was 0.019 and 0.045 across the forecasted range in Q_R, in the middle and lower segments, respectively. In the NewRE, surface pH_est always increased with Q_R, with ΔpH_est between 0.037 and 0.14.

FIGURE 3

Figure 3. Model pH_est as a function of atmospheric pCO₂ (μatm) and river discharge (m³ s^-1), where pH_est values are shown by the colored contours. Model results are shown for surface and bottom estuarine boxes, for the regions described in the Supplementary Material.

To summarize, the invasion of anthropogenic CO₂ resulted in moderate decreases in modeled pH_est, but this acidification was buffered when accompanied by elevated river discharge. The buffering effect of riverine inputs was greatest in surface waters and appears to be larger in the NewRE than in the NeuseRE. This difference between estuaries is particularly pronounced in the upper, oligohaline regions of both estuaries, where increases in Q_R were accompanied by rising pH in the NewRE but decreasing pH in the NeuseRE. It is important to note that this difference in pH change between surface and bottom waters may be overestimated by the box model, which does not consider mixing between these two layers. However, this model does provide an indication that the pH in both estuaries are sensitive to variations in fresh water delivery.

Inputs of TA from the watershed, in excess of DIC loads, may help buffer CO₂-driven acidification. Hence, variations in the ratio of TA:DIC have been used to indicate changes in the ability of water to buffer changes in pH (Hu and Cai, 2011). As shown in Figure 4, the TA:DIC ratio decreases exponentially with Q_R in both rivers. This suggests that both DIC and TA concentration in the river is primarily driven by weathering inputs, which are progressively diluted by surface water as Q_R increases. However, TA appears to decrease more rapidly than DIC as Q_R increases, causing the ratio of TA:DIC to decrease with rising Q_R, converging on a value of ∼0.5 at high discharge (Figure 4). At low Q_R, TA:DIC was relatively high in the New River, occasionally exceeding 1:1. TA:DIC was always below 1:1 in the Neuse River. Decadal trends in the TA:DIC ratio were stable in the NewRE, but were significantly negative (decreasing trend) in much of the upper NeuseRE (Supplementary Figure S3), corresponding with the decreasing pH trend in that estuary.

FIGURE 4

Figure 4. Relationship between river endmember TA:DIC and discharge for the New (red points) and Neuse (black points) rivers. Exponential fits are of the form f(x) = a^∗exp(b^∗x), where the coefficients are a = 0.74 and b = –0.0043 for the Neuse River, and a = 1.03 and b = –0.086 for the New River. Water quality data (TA, pH, Temperature, Salinity) were obtained from STORET for the North Carolina Dept. of Environmental Quality site J7850000 (Neuse River near Ft Barnwell), and from sites 21NC01WQ-P0600000 and 21NC01WQ-P1200000 in the upper New River. Only matching pairs of TA, pH, temp, and salinity were used, yielding a relatively small set of data from 1990 to 1993. DIC was calculated from measured TA, pH, Temperature, and Salinity with CO2SYS. Matching discharge data were obtained from the USGS gage at Gum Branch for the New River (2093000), and the USGS gage at Kinston for the Neuse River (2089500).

To better understand the effect of changing river end-member TA:DIC on buffering within the estuary, we calculated pH across the estuarine salinity gradient (with CO2SYS) for a variety of river end-member TA:DIC ratios that fall within the range of Figure 5A. This assessment is intended to indicate potential pH values for an idealized estuary, and does not include DIC/TA production or consumption by net metabolic processes or air-water CO₂ exchange. When salinity is greater than ∼20, carbonate chemistry is heavily controlled by mixing with well-buffered ocean water, thus, the decreasing riverine TA:DIC have limited effect on pH. However, when salinity falls below ∼10, the impact of variable riverine TA:DIC ratios becomes more apparent, causing large changes in calculated pH (Figure 5A). Furthermore, because the ratio of TA:DIC in the coastal ocean is greater than 1 (at least in this case), scenarios where riverine TA:DIC is less than 1 result in a point where TA:DIC passes through 1:1 zone, where buffering is locally minimized (Hu and Cai, 2013). To more precisely quantify this process, we calculate a buffer factor (β_DIC), which can be expressed in partial differential form as in Egleston et al. (2010):

${\mathrm{\beta }}_{\mathrm{D}\mathrm{I}\mathrm{C}}\mathrm{=}{\mathrm{\left(}{}^{\mathrm{\partial }\mathrm{I}\mathrm{n}\mathrm{\left[}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right]}\mathrm{/}}\mathrm{\partial }\mathrm{D}\mathrm{I}\mathrm{C}\mathrm{\right)}}^{\mathrm{-}\mathrm{1}}\mathrm{}\mathrm{\left(}\mathrm{6}\mathrm{\right)}$

FIGURE 5

Figure 5. Effect of varying river endmember TA:DIC ratios (1.1 = red, 1.0 = black, 0.9 = green, 0.8 = blue) on pH (A) and β_DIC (B, μmol kg^-1) within the estuary. The dashed lines show the effect of different H₂CO₃ dissociation constants on calculated pH_est [solid line: Millero, 2010 (M10), dash line: Cai and Wang, 1998 (CW98), and dotted line: Lueker et al., 2000 (L00)].

Expressed in this manner, increases in β_DIC indicate an increased buffering capacity, or resistance to changes in pH with DIC additions. As seen in Figure 5B, β_DIC exhibits a non-monotonic relationship with salinity (when TA:DIC is < 1), such that buffering is minimized at a salinity of ∼7 (riverine TA:DIC = 0.9) or ∼12 (TA:DIC = 0.8). Hence, when riverine TA:DIC is less than 1, estuarine pH is particularly sensitive to additional DIC inputs, which can enter via changes in net metabolism, anthropogenic CO₂ inputs, etc. This is the case for the NeuseRE, where TA:DIC appears to always be less than 1. That the greatest pH decreases in the NeuseRE occurred near the river mouth (-0.3 decade^-1, Figure 2A) further supports the role of poorly buffered river water on estuarine pH trends. TA:DIC was relatively high in the New River (Figure 4), causing downstream pH in the estuary to be less sensitive to Q_R. This concept is further supported by decadal trends in the TA:DIC ratio, which decreased at a significant rate across much of the upper NeuseRE, and was stable across the entire NewRE (Supplementary Figure S3). To summarize, it appears that TA delivery from the watershed (in excess of DIC) helped to buffer acidification stress in the NewRE, while decadal decreases in TA relative to DIC served to diminish the buffering capacity of the NeuseRE, allowing significant pH decreases to occur.

Summary of Acidification Drivers

A vast and tangled array of factors contribute to temporal variations in estuarine pH. Given this complexity, we determined in a statistically robust manner, the magnitude and direction of decadal-scale biogeochemical trends for two adjacent estuaries. Further, we identified which factors were the most important to the observed pH trends and determined if these factors differ between the two estuaries studied. While it is beyond the scope of this study to precisely quantify the many components affecting pH, we shed some light on the diverse drivers of estuarine acidification, a process that has gained interest in recent years.

Increases in river discharge [associated with ENSO activity (Supplementary Figure S2) during the Spring and Winter], coincided with a significant trend of acidification in the NeuseRE between 2005 and 2017. Acidification was greatest in bottom waters, concurrent with increased PPR in the surface, highlighting the importance of CO₂ inputs from net ecosystem respiration. Temporal trends in pH across the NewRE were small and not significant, despite moderate, but generally insignificant increases in phytoplankton biomass and productivity. Small, insignificant warming trends in both estuaries could not explain the observed bottom-water pH decrease. However, in the NewRE, this acidification ‘pressure’ was largely tempered by riverine loading of TA, which was higher relative to DIC. Inputs of TA from the Neuse River were smaller, however, leading to a weak buffer capacity that was more sensitive to autochthonous and allochthonous CO₂ inputs (Figure 5). While it is beyond the scope of this study to rigorously asses climate teleconnections like those between ENSO and annual-scale biogeochemical change, this topic may be of interest for future research, given the likely variations in ENSO with climate change.

Rising atmospheric CO₂ levels will act as a downward pressure on estuarine pH, but will be modified the net production or consumption of CO₂ by net ecosystem metabolism. Thus, in other stratified estuaries experiencing eutrophication, net respiratory inputs of CO₂ to bottom waters may interact with this anthropogenic CO₂ to geometrically lower pH. Meanwhile, concurrent change in river discharge (due to precipitation or human water consumption) have the potential to either enhance or ameliorate these acid inputs, although to an extent that varies according to differences in underlying lithology and hydrology. Hence, all estuaries experiencing future changes in river discharge may experience similar pH trajectories due to the interaction between atmospheric CO₂ inputs and riverine hydrogeochemistry. Local decision makers and stakeholders should consider these factors when considering management strategies aimed at coping with future estuarine acidification.

Data Availability

The DCERP (https://dcerp.rti.org/) and ModMon (http://paerllab.web.unc.edu/projects/modmon/products-and-applications/) datasets are published on the data sharing repository Figshare (https://figshare.com/articles/WQ_Datasets_-_Decadal_Acidification_Trends/7643195).

Author Contributions

BVD and HW contributed equally to the conception and design of the study.

Funding

This research was funded by SERDP-DCERP project RC-2245, the North Carolina Department of Environmental Quality (ModMon Program), the Lower Neuse Basin Association, NC Sea Grant, and the UNC Water Resources Research Institute.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Hans Paerl, Betsy Abare, Jeremey Braddy, Karen Rossignol, and Randy Sloup of the University of North Carolina at Chapel Hill, Institute of Marine Sciences, for compiling and maintaining ModMon and DCERP databases. Xinping Hu offered valuable guidance during the preparation of this manuscript. We also thank the organizers of the 2017 Coastal and Estuarine Research Federation conference, at which the concept for this manuscript was conceived.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2019.00136/full#supplementary-material

References

Aufdenkampe, A. K., Mayorga, E., Raymond, P. A., Melack, J. M., Doney, S. C., Alin, S. R., et al. (2011). Riverine coupling of biogeochemical cycles between land, oceans, and atmosphere. Front. Ecol. Environ. 9:53–60. doi: 10.1890/100014

CrossRef Full Text | Google Scholar

Bates, N., Astor, Y., Church, M., Currie, K., Dore, J., and Gonaález-Dávila, M. (2014). A time-series view of changing ocean chemistry due to ocean uptake of anthropogenic CO2 and ocean acidification. Oceanography 27, 126–141. doi: 10.5670/oceanog.2014.16

CrossRef Full Text | Google Scholar

Baumann, H., and Smith, E. M. (2017). Quantifying metabolically driven pH and oxygen fluctuations in US nearshore habitats at diel to interannual time scales. Estuaries Coasts 41, 1102–1117. doi: 10.1007/s12237-017-0321-3

CrossRef Full Text | Google Scholar

Boyer, J. N., Christian, R. R., and Stanleyl, D. W. (1993). Patterns of phytoplankton primary productivity in the Neuse River estuary, North Carolina, USA. Mar. Ecol. Prog. Ser. 97, 287–297. doi: 10.3354/meps097287

CrossRef Full Text | Google Scholar

Cai, W. J., Hu, X., Huang, W. J., Murrell, M. C., Lehrter, J. C., and Lohrenz, S. E. (2011). Acidification of subsurface coastal waters enhanced by eutrophication. Nat. Geosci. 4, 766–770. doi: 10.1038/ngeo1297

CrossRef Full Text | Google Scholar

Cai, W. J., Huang, W. J., Luther, G. W., Pierrot, D., Li, M., and Testa, J. (2017). Redox reactions and weak buffering capacity lead to acidification in the Chesapeake Bay. Nat. Commun. 8:369. doi: 10.1038/s41467-017-00417-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, W.-J., and Wang, Y. (1998). The chemistry, fluxes, and sources of carbon dioxide in the estuarine waters of the Satilla and Altamaha Rivers, Georgia. Limnol. Oceanogr. 43, 657–668. doi: 10.4319/lo.1998.43.4.0657

CrossRef Full Text | Google Scholar

Couldrey, M. P., Oliver, K. I. C., Yool, A., Halloran, P. R., and Achterberg, E. P. (2016). On which timescales do gas transfer velocities control North Atlantic CO2 flux variability? Global biogeochem. Cycles 30, 787–802. doi: 10.1002/2015GB005267

CrossRef Full Text | Google Scholar

Crosswell, J. R., Anderson, I. C., Stanhope, J. W., Dam, B. V., Brush, M. J., and Ensign, S. (2017). Carbon budget of a shallow, lagoonal estuary: transfor-sink dynamics along the river-estuary-ocean continuum. Limnol. Oceanogr. 62, S29–S45. doi: 10.1002/lno.10631

CrossRef Full Text | Google Scholar

Crosswell, J. R., Paerl, H. W., Wetz, M. S., and Hales, B. (2012). Air-water CO2 fluxes in the microtidal Neuse River Estuary, North Carolina. J. Geophys. Res. 117:12. doi: 10.1029/2012JC007925

CrossRef Full Text | Google Scholar

Cyronak, T., Andersson, A. J., D’Angelo, S., Bresnahan, P., Davidson, C., and Griffin, A. (2018). Short-term spatial and temporal carbonate chemistry variability in two contrasting seagrass meadows: implications for pH buffering capacities. Estuaries Coasts 41, 1282–1296. doi: 10.1007/s12237-017-0356-5

CrossRef Full Text | Google Scholar

Doney, S. C., Fabry, V. J., Feely, R. A., and Kleypas, J. A. (2009). Ocean acidification: the other CO2 problem. Ann. Rev. Mar. Sci. 1, 169–192. doi: 10.1146/annurev.marine.010908.163834

CrossRef Full Text | Google Scholar

Doney, S. C., and Schimel, D. S. (2007). Carbon and climate system coupling on timescales from the precambrian to the anthropocene. Annu. Rev. Environ. Resour. 32, 31–66. doi: 10.1146/annurev.energy.32.041706.124700

CrossRef Full Text | Google Scholar

Egleston, E. S., Sabine, C. L., and Morel, F. M. M. (2010). Revelle revisited: buffer factors that quantify the response of ocean chemistry to changes in DIC and alkalinity. Global Biogeochem. Cycles 24, 1–9. doi: 10.1029/2008GB003407

CrossRef Full Text | Google Scholar

Fear, J. M., Paerl, H. W., and Braddy, J. S. (2007). Importance of submarine groundwater discharge as a source of nutrients for the Neuse River Estuary, North Carolina. Estuaries Coasts 30, 1027–1033. doi: 10.1007/BF02841393

CrossRef Full Text | Google Scholar

Feely, R. A., Alin, S. R., Newton, J., Sabine, C. L., Warner, M., Devol, A., et al. (2010). The combined effects of ocean acidification, mixing, and respiration on pH and carbonate saturation in an urbanized estuary. Estuar. Coast. Shelf Sci. 88, 442–449. doi: 10.1016/j.ecss.2010.05.004

CrossRef Full Text | Google Scholar

Hall, N. S., Paerl, H. W., Peierls, B. L., Whipple, A. C., and Rossignol, K. L. (2012). Effects of climatic variability on phytoplankton community structure and bloom development in the eutrophic, microtidal, New River Estuary, North Carolina, USA. Estuar. Coast. Shelf Sci. 117, 70–82. doi: 10.1016/j.ecss.2012.10.004

CrossRef Full Text | Google Scholar

Herrmann, M., Najjar, R. G., Kemp, W. M., Alexander, R. B., Boyer, E. W., Cai, W.-J., et al. (2015). Net ecosystem production and organic carbon balance of U.S. East Coast estuaries: a synthesis approach. Global Biogeochem. Cycles 29, 96–111. doi: 10.1002/2013GB004736

CrossRef Full Text | Google Scholar

Hu, X., and Cai, W. J. (2011). An assessment of ocean margin anaerobic processes on oceanic alkalinity budget. Glob. Biogeochem. Cycles 25, 1–11. doi: 10.1029/2010GB003859

CrossRef Full Text | Google Scholar

Hu, X., and Cai, W. J. (2013). Estuarine acidification and minimum buffer zone. A conceptual study. Geophys. Res. Lett. 40, 5176–5181. doi: 10.1002/grl.51000

CrossRef Full Text | Google Scholar

Hu, X., Pollack, J. B., McCutcheon, M. R., Montagna, P. A., and Ouyang, Z. (2015). Long-term alkalinity decrease and acidification of estuaries in northwestern Gulf of Mexico. Environ. Sci. Technol. 49, 3401–3409. doi: 10.1021/es505945p

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeffrey, L. C., Maher, D. T., Santos, I. R., McMahon, A., and Tait, D. R. (2016). Groundwater, acid and carbon dioxide dynamics along a coastal wetland, lake and estuary continuum. Estuaries Coasts 39, 1325–1344. doi: 10.1007/s12237-016-0099-8

CrossRef Full Text | Google Scholar

Keeling, R. F., and Keeling Charles, D. (2017). Atmospheric Monthly In Situ CO2 Data - MaunaLoa Observatory, Hawaii. In Scripps CO2 Program Data. UC San Diego Library Digital Collections. Available at: http://doi.org/10.6075/J08W3BHW

Le Quéré, C., Andrew, R. M., Canadell, J. G., Sitch, S., Korsbakken, J. I., and Peters, G. P. (2016). Global Carbon Budget 2016. Earth Syst. Sci. Data Discuss. 8, 605–649. doi: 10.5194/essd-8-605-2016

CrossRef Full Text | Google Scholar

Lewis, E., and Wallace, D. (1998). Program Developed for CO2 System Calculations. Oak Ridge, TN: Carbon Dioxide Information Analysis Center. doi: 10.2172/639712

CrossRef Full Text | Google Scholar

Lueker, T. J., Dickson, A. G., and Keeling, C. D. (2000). Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity and equations K1 and K2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium. Mar. Chem. 70, 105–119. doi: 10.1016/S0304-4203(00)00022-0

CrossRef Full Text | Google Scholar

Mallin, M., Paerl, H. W., Rudek, J., and Bates, P. (1993). Regulation of estuarine primary production by watershed rainfall and river flow. Mar. Ecol. Prog. Ser. 93, 199–203. doi: 10.3354/meps093199

CrossRef Full Text | Google Scholar

Mallin, M. A., McIver, M. R., Wells, H. A., Parsons, D. C., and Johnson, V. L. (2005). Reversal of eutrophication following sewage treatment upgrades in the New River Estuary. North Carolina. Estuaries 28, 750–760. doi: 10.1007/BF02732912

CrossRef Full Text | Google Scholar

Millero, F. J. (2010). Carbonate constants for estuarine waters. Mar. Freshw. Res. 61, 139–142. doi: 10.1071/MF09254

CrossRef Full Text | Google Scholar

Mote, P. W., Parson, E., Hamlet, A. F., Keeton, W. S., Lettenmaier, D., and Mantua, N. (2008). Preparing for climatic change: the water, salmon, and forests of the Pacific Northwest. Clim. Change 61, 45–88. doi: 10.1023/A:1026302914358

CrossRef Full Text | Google Scholar

Müller, J. D., Schneider, B., and Rehder, G. (2016). Long-term alkalinity trends in the Baltic Sea and their implications for CO2-induced acidification. Limnol. Oceanogr. 61, 1984–2002. doi: 10.1002/lno.10349

CrossRef Full Text | Google Scholar

Nixon, S. W., Oczkowski, A. J., Pilson, M. E. Q., Fields, L., Oviatt, C. A., Hunt, C. W., et al. (2015). On the response of pH to inorganic nutrient enrichment in well-mixed coastal marine waters. Estuaries Coasts 38, 232–241. doi: 10.1007/s12237-014-9805-6

CrossRef Full Text | Google Scholar

Null, K. A., Corbett, D. R., DeMaster, D. J., Burkholder, J. M., Thomas, C. J., Reed, R. E., et al. (2011). Porewater advection of ammonium into the Neuse River Estuary, North Carolina, USA. Estuar. Coast. Shelf Sci. 95, 314–325. doi: 10.1016/j.ecss.2011.09.016

CrossRef Full Text | Google Scholar

Pacella, S. R., Brown, C. A., Waldbusser, G. G., Labiosa, R. G., and Hales, B. (2018). Seagrass habitat metabolism increases short-term extremes and long-term offset of CO2 under future ocean acidification. PNAS 115, 1–6. doi: 10.23719/1407616

PubMed Abstract | CrossRef Full Text | Google Scholar

Paerl, H. W., Pinckney, J., Fear, J., and Peierls, B. L. (1998). Ecosystem responses to internal and watershed organic matter loading: consequences for hypoxia in the eutrophying Neuse River Estuary, North Carolina, USA. Mar. Ecol. Prog. Ser. 166, 17–25. doi: 10.3354/meps166017

CrossRef Full Text | Google Scholar

Peierls, B. L., Hall, N. S., and Paerl, H. W. (2012). Non-monotonic responses of phytoplankton biomass accumulation to hydrologic variability: a comparison of two coastal plain north carolina estuaries. Estuaries Coasts 35, 1376–1392. doi: 10.1007/s12237-012-9547-2

CrossRef Full Text | Google Scholar

Phlips, E. J., Badylak, S., Hart, J., Haunert, D., Lockwood, J., and O’Donnell, K. (2011). Climatic influences on autochthonous and allochthonous phytoplankton blooms in a subtropical estuary, St. Lucie Estuary, Florida, USA. Estuaries Coasts 35, 335–352. doi: 10.1007/s12237-011-9442-2

CrossRef Full Text | Google Scholar

Reimer, J. J., Wang, H., Vargas, R., and Cai, W.-J. (2017). Multidecadal fCO2 increase along the United States Southeast Coastal Margin. J. Geophys. Res. Ocean 122, 10061–10072. doi: 10.1002/2017JC013170

CrossRef Full Text | Google Scholar

Siam, M. S., and Eltahir, E. A. B. (2017). Climate change enhances interannual variability of the Nile river flow. Nat. Clim. Chang. 7, 350–354. doi: 10.1038/nclimate3273

CrossRef Full Text | Google Scholar

Takahashi, T., Sutherland, S. C., Chipman, D. W., Goddard, J. G., and Ho, C. (2014). Climatological distributions of pH, pCO2, total CO2, alkalinity, and CaCO3saturation in the global surface ocean, and temporal changes at selected locations. Mar. Chem. 164, 95–125. doi: 10.1016/j.marchem.2014.06.004

CrossRef Full Text | Google Scholar

Unsworth, R. K. F., Collier, C. J., Henderson, G. M., and McKenzie, L. J. (2012). Tropical seagrass meadows modify seawater carbon chemistry: implications for coral reefs impacted by ocean acidification. Environ. Res. Lett. 7:024026. doi: 10.1088/1748-9326/7/2/024026

CrossRef Full Text | Google Scholar

Van Dam, B. R., Crosswell, J. R., Anderson, I. C., and Paerl, H. W. (2018a). Watershed-Scale Drivers of Air-Water CO2Exchanges in Two Lagoonal North Carolina (USA) Estuaries. J. Geophys. Res. Biogeosci. 123, 271–287. doi: 10.1002/2017JG004243

CrossRef Full Text | Google Scholar

Van Dam, B. R., Crosswell, J. R., and Paerl, H. W. (2018b). Flood-driven CO2 emissions from adjacent North Carolina estuaries during Hurricane Joaquin (2015). Mar. Chem. 207, 1–12. doi: 10.1016/j.marchem.2018.10.001

CrossRef Full Text

Waldbusser, G. G., and Salisbury, J. E. (2014). Ocean acidification in the coastal zone from an organism’s perspective: multiple system parameters, frequency domains, and habitats. Ann. Rev. Mar. Sci. 6, 221–247. doi: 10.1146/annurev-marine-121211-172238

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Hu, X., Cai, W. J., and Sterba-Boatwright, B. (2017). Decadal fCO2 trends in global ocean margins and adjacent boundary current-influenced areas. Geophys. Res. Lett. 44, 8962–8970. doi: 10.1002/2017GL074724

CrossRef Full Text | Google Scholar

Wang, H., Hu, X., and Sterba-Boatwright, B. (2016). A new statistical approach for interpreting oceanic fCO2 data. Mar. Chem. 183, 41–49. doi: 10.1016/j.marchem.2016.05.007

CrossRef Full Text | Google Scholar

Wang, Z. A., Kroeger, K. D., Ganju, N. K., Gonneea, M. E., and Chu, S. N. (2016). Intertidal salt marshes as an important source of inorganic carbon to the coastal ocean. Limnol. Oceanogr. 61, 1916–1931. doi: 10.1002/lno.10347

CrossRef Full Text | Google Scholar