Contrasting trends in water quality between adjacent ocean‐ and river‐dominated estuaries: Evidence for marsh porewaters as a source of nutrient enrichment?

Degradation of estuarine water quality during the Anthropocene has largely resulted from discharges of nutrients leading to eutrophication. Recently, upstream management practices have led to comparatively reduced nutrient input into estuaries. Concurrently, climate cycles and impacts associated with anthropogenic climate warming can affect the long‐term conditions observed within estuaries. Using long‐term monitoring data from adjacent southeastern U.S. estuaries, we show that decadal‐scale trends in nutrient concentrations and phytoplankton standing stock differ between the two connected systems. These contrasting trends appear to result from differences in oceanic influence, the extent of adjacent vegetated marsh, watershed size, and upstream degradation. In the minimally impacted, ocean‐dominated North Inlet estuary, we document increasing ammonium and chlorophyll a (Chl a), while in the adjacent, river‐dominated Winyah Bay, ammonium, and Chl a concentrations are more variable but do not appear to have increased over the same time period. Surprisingly, total nitrogen exhibits the opposite pattern: temporal stability in North Inlet but increasing in Winyah Bay. We hypothesize that sea level rise associated with climate change has driven a complex set of interactions between salt marsh porewaters and tidal pumping, leading to the spillover of nutrients from salt marshes into tidal creeks in North Inlet. In Winyah Bay, this mechanism is less evident as a driver of ammonium concentrations, likely due to the outsized effect of watershed nutrient input and the narrow fringing marsh platform. The degree to which this mechanism operates in other estuaries, which vary in tidal range, the extent of vegetated marsh, watershed size, and degree of anthropogenic degradation warrants further study.

However, alongside anthropogenic eutrophication, climate also plays a role in nutrient cycling and the long-term conditions observed within estuarine ecosystems (Rabalais et al. 2009).For example, in the Chesapeake Bay, precipitation and tropical cyclones lead to variability in primary production which is overlaid on longer-term increases in chlorophyll a (Chl a) driven by eutrophication (Harding et al. 2016).Similarly, in a tributary of the Chesapeake, climate warming is predicted to increase hypoxia despite concurrent reductions in nutrient concentrations (Testa et al. 2022).The mixed effects of climate and nutrient loading on phytoplankton dynamics result in a complex signal making it difficult to characterize long-term trends in water quality.Nonetheless, quantification of these trends is invaluable for understanding the relative importance of climate effects and pollution-reduction strategies on estuarine water quality.In addition, comparing trends in nutrients and phytoplankton between different estuaries with variable degrees of anthropogenic impact would provide insights into the range of responses for coastal ecosystems facing similar regional climate forcing.
Despite rapid suburban development and a continued increase in human population growth along much of the U.S. coastline (Freeman et al. 2019), and the southeastern U.S. in particular (Chi and Ho 2018), pockets of coastal habitat remain undeveloped and their adjacent estuaries can appear, at least superficially, minimally impacted.Nevertheless, highfrequency, long-term sampling of environmental parameters within these minimally impacted ecosystems can reveal unexpected conditions with respect to water quality.For example, even with low levels of development in the surrounding watershed, tidal creeks within the Ashepoo-Combahee-Edisto Basin in South Carolina exhibited elevated nutrient and Chl a concentrations relative to other similar estuaries (Keppler et al. 2015).Similarly, data collected on a suite of environmental conditions across estuarine systems nationwide revealed that atmospheric and estuarine heatwaves commonly co-occur and can be accompanied by periods of low dissolved oxygen and/or low pH (Tassone et al. 2021), representing a possible triple threat to even minimally impacted estuarine systems and the organisms which rely on them.Given the ongoing and predicted future stressors to coastal ecosystems, rigorous long-term monitoring programs provide a basis to understand the current state of estuaries, the effects of previous disturbance events or management decisions, and the changes driven by upcoming human impacts.However, quantifying the trends and trajectories in these data has proven difficult due to the multiple scales of variability, including tidal, meteorological, annual, and climatic cycles which range from hours to decades in time scale, and which are overlaid on changes in nutrient loading and global climate warming.
The National Estuarine Research Reserve (NERR) System is a state-federal partnership which operates a network of 30 coastal reserves across the United States.Each NERR conducts long-term monitoring of environmental conditions through the System-Wide Monitoring Program (SWMP), which is designed to track short-term variability and longterm change in estuarine water quality and weather.SWMP data collection efforts include water temperature, salinity, pH, dissolved oxygen, and concentrations of nutrients and Chl a, among other parameters (Buskey et al. 2015;Kennish 2019).A previous analysis of SWMP data collected within the North Inlet-Winyah Bay (NI-WB) NERR, located along the South Carolina coast, during the period of 1993-2001 demonstrated that levels of dissolved oxygen, N-and P-based nutrients, and Chl a were within the range of values previously reported from other southeastern U.S. estuaries (Buzzelli et al. 2004).However, the river-dominated and ocean-dominated subsystems which comprise the Reserve exhibited starkly different levels of Chl a, salinity, and some nutrient constituents (e.g., dissolved organic carbon, nitrate, and orthophosphate).Furthermore, Buzzelli et al. (2004) did not consider trends over time within or between the two estuaries.Thus, to date, we lack an understanding of the potential differences in long-term trends in water quality between these adjacent but contrasting ecosystems: an ocean-dominated and river-dominated estuary.In addition, since 2001, sea level rise has accelerated along the coastlines of the southeastern U.S., with continued increases in projected sea level heights of $ 0.4 m by 2050 and up to $ 2 m by 2150 (Sweet et al. 2022).The NI-WB system represents an opportunity to investigate changes in water quality conditions through time within two connected but dissimilar estuaries undergoing rapid sea level rise.
Here, we utilized nearly two decades of SWMP water quality data collected at the NI-WB NERR to quantify the long-term trends and periodic variability in nutrients and phytoplankton between two connected estuaries which differ in oceanic and watershed influence as well as upstream anthropogenic alteration.We hypothesize that sea level rise associated with climate change has driven decadal-scale increases in inorganic nutrient concentrations within the oceandominated North Inlet estuary.We suggest that this is primarily due to complex interactions between salt marsh porewaters and tidal pumping associated with increasing inundation time (Krask et al. 2022).Conversely, we propose that riverine input from the large Winyah Bay watershed exerts overwhelming control over nutrient concentrations in the Winyah Bay estuary, such that any changes resulting from sea level rise are comparatively minimal.Our results demonstrate how water quality within ocean-and river-dominated estuaries can respond differently to rising sea level based on the relative importance of several drivers of nutrient cycling.This may provide insight into patterns emerging in similar estuarine environments at temperate latitudes.

Study area and water quality monitoring data collection
The North Inlet-Winyah Bay National Estuarine Research Reserve (NI-WB NERR) is located in Georgetown, South Carolina, USA and encompasses two connected, yet very different, estuarine systems.North Inlet estuary is a relatively small (33 km 2 ), bar-built, ocean-dominated, high salinity estuary (Fig. 1).Much of North Inlet estuary consists of Spartina alterniflora-dominated salt marsh, as well as well-mixed tidal creeks that drain a small (38 km 2 ), largely forested watershed.Freshwater input into the North Inlet estuary is limited to run-off from the surrounding watershed and infrequent pulses from adjacent Winyah Bay through connected tidal creeks during periods of high river discharge (Traynum and Styles 2008).The North Inlet watershed contains minimal development, all of which primarily occurred decades before the time period of our current study.Thus, any recent trends in environmental parameters within the estuary are likely due to larger-scale processes rather than direct anthropogenic impacts on the watershed.In contrast, Winyah Bay is a large (65 km 2 ), coastal plain, river-dominated system with variable salinity and freshwater inflow from five rivers.The watershed for Winyah Bay drains $ 47,000 km 2 of predominantly agricultural land (Allen et al. 2014; Fig. 1).Thus, our two adjacent study sites represent substantially different types of estuarine ecosystems.Within each of these two estuaries, we focused on a single, representative sampling location: Clambank Landing (CB) in the center of North Inlet, and Thousand Acre (TA) located at the mouth of a tidal creek immediately adjacent to the mainstem of Winyah Bay.Our analysis focuses on key parameters used to characterize water quality status within estuaries (ammonium, nitrate, total nitrogen, and Chl a concentrations) as well as the ratio of dissolved inorganic nitrogen to dissolved inorganic phosphorus, DIN : DIP.Ammonium (NH 4 ) is a key inorganic nutrient constituent assimilated by phytoplankton in primary production (Dugdale et al. 2007;Domingues et al. 2011).Nitrate (NO 3 ), another inorganic nitrogen species, concentrations are strongly affected by land use patterns within a watershed (Justi c et al. 2003) and are generally positively related to human development and agriculture (Basnyat et al. 1999).However, these inorganic forms of nitrogen represent only a small fraction of the overall pool in estuarine systems.Total nitrogen (TN) accounts for the inorganic nitrogen fraction as well as the larger pool of organic compounds, some of which may be bioavailable to microorganisms or can be recycled into more readily usable forms of nitrogen (i.e., NH 4 ) upon organic matter decomposition and remineralization (Twomey et al. 2005).Algal biomass (measured as Chl a) provides a measure of ecosystem productivity and an index of autotrophic carbon production (Cloern 2001).We define DIN as the sum of ammonium, nitrate, and nitrite concentrations, while DIP is the concentration of orthophosphate (PO 4 ).
Water samples were collected on a 20-d interval during the period from 2002 to 2019 at CB (33.333883 N, 79.192923 W) and TA (33.298252 N, 79.256473 W).This sampling interval was selected to avoid systematic bias associated with the spring-neap tidal periodicity observed in semi-diurnal systems while balancing logistic feasibility for a long-term sampling program (Kjerfve and McKellar 1980;Kjerfve and Wolaver 1988).During a given sampling event at each station, automated water samplers (ISCO 6712) were pre-programmed to collect 950 mL samples from a depth of 0.5 m over the course of two complete tidal cycles (24 h and 48 min), beginning at the first low tide of the day and occurring subsequently every 2 h and 4 min (13 samples total).
Over the course of each sampling event, samples were kept in the ISCOs (due to limited site accessibility), which were preloaded with frozen, gallon-sized water jugs at the core of the ISCO base to keep samples cold.Following the collection of the final sample of the event, the ISCOs were retrieved, and samples were immediately stored at 4 C as they were filtered (Whatman GF/F; nominal pore size of 0.7 μm) in batches over the subsequent 4-6 h.Filtered samples were used in analyses for all nutrient parameters except for TN, for which an aliquot of whole (unfiltered) water was taken before sample filtration.Samples remained stored at 4 C and were analyzed for NH 4 , NO 2 , NO x (= NO 2 + NO 3 ), and PO 4 the following day (within 48 h of collection), while the unfiltered sample was analyzed for TN within 28 d (see Supporting Information Appendix S2 for analytical laboratory procedures).Despite the addition of ice to the ISCOs, we acknowledge the risk of potential changes resulting from ongoing nutrient cycling processes occurring between the times of sample collection, filtration, and analysis.However, we are confident that any interim changes are negligible because samples collected over the $ 25 h tidal cycle tend to show a clear high-low tidal signal.Specifically, NO x and NH 4 concentrations between recurring tidal stages within a sampling event are typically highly similar despite their different functional hold times and we see no evidence of a systematic directional trend within sampling events throughout the time series (Supporting Information Fig. S1).This suggests that any interim changes in nutrient concentrations due to sample hold times are negligible compared with variability associated with other environmental drivers.
In addition to the standard quality control methods utilized by the NERR Centralized Data Management Office (CDMO: cdmo.baruch.sc.edu/data/qaqc.cfm),we conducted additional QA/QC and data pooling before fitting models in our analyses.Specifically, we calculated the means of the consecutive high tide (n = 2) and low tide (n = 3) samples at a given site for each parameter for a particular sampling event (every 20 d).From those event-level mean values, any data point that was > 5 standard deviations from the overall, longterm mean was excluded from further analysis.Samples that were below the defined minimum detection limit (MDL) were reported as the MDL value (Supporting Information Appendix S2).These data points were retained in the data set so as not to bias trends toward higher values.

Analytical methods: Time series trends
We first aimed to understand the relative proportions of the available forms of nitrogen at both sampling stations and how they may have changed over time.For the periods 2002-2004 and 2017-2019, we used measured values of NO 2 + NO 3 , NH 4 , total dissolved nitrogen, and TN to calculate monthly mean concentrations of NO 2 + NO 3 , NH 4 , dissolved organic nitrogen (DON) and particulate nitrogen (PN).In sum, these four constituents are equal to the TN pool.
To illustrate and analyze trends in our water quality time series data (mean high and low tide samples collected at a 20-d interval), we utilize generalized additive models (GAMs), an analytical method emerging for use within estuarine science (Testa et al. 2018;Murphy et al. 2019;Beck et al. 2022).These models can effectively characterize non-linear time series trends and account for predictors operating at different scales (i.e., intra-annual/seasonal patterns).Generally, a GAM is a regression-based statistical model that allows for nonlinear fits based on smoothed functions that take their shape from the underlying data (Simpson 2018).GAMs produce comparable results to other modeling approaches often used for water quality applications (Beck and Murphy 2017).In our case, we used GAMs for non-linear time series trend estimation as well as identifying periods of change in the trend using the first derivative.For each of the five focal parameters measured at both sampling stations-NH 4 ; NO 3 (calculated as NO x -NO 2 ); TN; DIN : DIP; and Chl a-we fit a GAM with terms s("day") of the time series and s("month") of year, the latter of which can account for intra-annual or seasonal patterns.For all models, the basis complexity (k) was set to 15 for s(day) and 12, for s(month).In addition, we set the smooth term for Month to be cyclical such that the end-points of the spline are forced to be equal (by setting bs = "cc" in our call to the gam function).
From each fitted GAM, we obtain estimates of the smoothed effects for the time series trend and intra-annual trend (based on their respective estimated degrees of freedom, a measurement of "wiggliness"), fitted values, residuals, goodness-of-fit statistics (adjusted-R 2 and % deviance explained), and associated F statistics and p values.Models were fit via restricted maximum likelihood.For each GAM, the response variable was a nutrient concentration in mg L À1 or Chl a concentration in μg L À1 .To characterize points in time at which the slope of the time series trend changes, we utilized the first derivative technique of Simpson (2018).This method uses the simultaneous confidence intervals around the estimate of the first derivative of the partial effect of day of the time series; periods of statistically significant change in the time series trend are identified as those time periods where these first derivative confidence intervals do not include zero.

Analytical methods: Time series periodic components
To investigate trends in periodicity in water quality parameters within NI-WB over the period 2002-2019, we conducted wavelet analyses.Spectral analysis using wavelets can extract periodic components from non-stationary time series (Cazelles et al. 2008), and a key benefit of this approach for environmental time series is the ability to observe how periodicities change over the course of the time series (Winder and Cloern 2010;Rasmuson and Shanks 2020).We used the continuous wavelet transform with a Morlet wavelet (wavenumber, ω 0 = 6; Torrence and Compo 1998) to investigate periodicities within the power spectra of each parameter/ tide height/sampling site time series.Any missing sampling events were interpolated from the neighboring data.In addition, we calculated time-averaged ("global") wavelet power spectra, along with their 95% significance levels (based on a red noise background spectrum), to visualize the dominant periodicities across each time series.

Analytical methods: Drivers of ammonium across estuary types
Finally, to compare the relative importance of different potential drivers of nutrient cycling within our two focal estuaries, we conducted random forest analyses on annual mean values of NH 4 concentration at both high and low tide.Random forest is a nonparametric classification method that uses random subsets of the dataset to quantify parameter importance, and thus the ability of each independent variable to predict the response variable (Breiman 2001).Based on the trees produced by random forest, an importance value is assigned to each predictor variable, and we normalized importance values to sum to 1 because relative, not absolute, importance values are informative.Typically, random forest is utilized on datasets with many potential predictors.However, considering that this framework provides an intuitive way to compare the relative strength of predictor variables, we have adopted it here.Thus, given our relatively small dataset (18 annual observations and 5 potential drivers) and to account for the stochastic nature of random forest resampling, we chose to run 1000 random forest simulations for each site/ tidal stage combination.From these, we then calculated the mean and standard error of the importance values associated with each predictor variable.In our case, predictors with large importance values can be considered as key correlates of NH 4 concentration (our response variable), with separate analyses conducted for high and low tide data from each of the TA and CB monitoring stations.We first calculated station-specific mean NH 4 concentrations for each year of our time series, 2002-2019 (Supporting Information Fig. S2).We then collated data describing five potential drivers of water quality in our focal system: (1) river flow in the Pee Dee River, which is the major tributary into Winyah Bay; (2) sea level height from a NOAA tide gauge located at Oyster Landing, within the North Inlet estuary; (3) precipitation at the NI-WB NERR SWMP meteorological station at Oyster Landing, in the North Inlet estuary; (4) the multivariate El Niño Southern Oscillation Index (MEI), which is based on multiple variables that holistically describe conditions in the tropical Pacific during ENSO events; and (5) the North Atlantic Oscillation (NAO) index, which characterizes differences in surface sea level pressure between the subtropical and subpolar Atlantic.Specifically, Pee Dee River discharge data (ft 3 s À1 ) were retrieved from the United States Geological Survey monitoring station at Pee Dee, South Carolina (USGS 02131000), which is located upstream from Winyah Bay.Discharge measurements taken at 15-or 30-min intervals (variable across the study period) were used to calculate a mean annual river flow for use in our random forest analysis (Supporting Information Fig. S3).Monthly sea level height data were retrieved from the National Oceanic and Atmospheric Administration tide gauge at Oyster Landing (NOAA station 8662245) then aggregated to calculate a mean annual sea level height (Supporting Information Fig. S4).Precipitation data were retrieved from the NERRS CDMO, similar to the nutrient data described above and aggregated to calculate annual mean precipitation (Supporting Information Fig. S4).MEI data were retrieved from NOAA's National Center for Atmospheric Research: https://psl.noaa.gov/data/correlation/meiv2.data.MEI is initially computed for 12 bimonthly sliding windows per year; we aggregated these to calculate the annual mean for each year of our study period (Supporting Information Fig. S3).Finally, NAO index data were retrieved from NOAA's National Center for Environmental Information: https://www.ncdc.noaa.gov/teleconnections/nao/(Supporting Information Fig. S3).
All analyses were conducted in R version 4.0.5.GAMs were constructed with the mgcv package (Wood 2021).Wavelet and power spectra analyses were conducted with the biwavelet package (Gouhier et al. 2019).Random forest analysis was based on the protocols of Harper et al. (2011) and conducted using the randomForest package (Breiman et al. 2015).

Time series trends
The monthly mean averages of N constituents for the periods 2002-2004 and 2017-2019 demonstrate that overall concentrations of all forms were generally higher at TA relative to CB throughout the year at the beginning and end of the time series (Fig. 2).For both time periods, NH 4 is the dominant inorganic form of nitrogen present at CB, while the proportion and overall concentration of NO x in the inorganic N fraction is much higher at TA.However, the inorganic constituents generally comprise a much smaller proportion of the TN pool relative to the DON and PN forms at both sites for both time periods.Average monthly TN concentrations tend to be higher for a given month at TA during the 2017-2019 time frame compared with 2002-2004, while TN concentrations at CB are consistent across the two periods.Within the TN pool, monthly average NH 4 concentrations appear to have increased at CB in the latter period while the relative proportion of DON has decreased.Finally, NO x concentrations are generally higher at TA from 2017 to 2019 compared with 2002-2004.
Overall, our GAMs provided a good fit to nutrient and Chl a time series data collected at CB, with the exception of nitrate, NO 3 .Models effectively exhibited both the overall trends for 2002-2019 and the expected intra-annual, seasonal patterns (Figs.3-5; Supporting Information Fig. S5; Table 1).Excluding NO 3 , models fit to data from CB explained an average of 35% of deviance for nutrient parameters across both high and low tide heights (range: 22-48%) and explained > 60% of deviance for Chl a data collected at low tide (Table 1).GAMs fit poorly to NO 3 data from CB, only explaining $ 5% of deviance at either high or low tide (Table 1).Models describing nutrient time series data from TA fit less well for nutrient parameters (excluding NO 3 ), explaining on average 28% of deviance and as low as 15% of deviance for TN and DIN : DIP at high and low tide, respectively (Table 1).However, unlike NO 3 at CB, GAMs explained a substantial amount of deviance for NO 3 data from TA (56% and 42% deviance explained for high and low tide, respectively).With four exceptions (NO 3 at both tides, Chl a at high tide, and TN at low tide), the model fits were worse (in terms of percent deviance explained) for TA compared with CB (Table 1).
Our models also demonstrate the differing trends in most nutrients and Chl a time series between CB and TA.Generally, TA exhibited higher nutrient and Chl a concentrations, as well as more variability with the little directional trend, while the North Inlet station at CB had substantially lower variability and in some cases increased over time (Figs.3-5; Supporting Information Fig. S5).For example, NH 4 at CB increased from 2002 to 2019, with a statistically significant jump in the slope of the trend during 2011 (Fig. 3).Similarly, as a result of the increase in NH 4 but relatively stable PO 4 over time (Supporting Information Fig. S5), the DIN : DIP ratio also increased through time at CB (Supporting Information Fig. S5).Chl a concentration also significantly increased at CB during our study period, but the increase in the slope of the trend occurred with a 4-yr time lag compared to the increase in NH 4 (Fig. 4).Conversely, TA did not display statistically significant increases in the slope of the trend in NH 4 , DIN : DIP, or Chl a between 2002 and 2019 (Supporting Information Fig. S6).
A key exception, however, is that of NO 3 which showed no change point at CB but exhibited a strong increase over time at TA (Supporting Information Fig. S7).Despite the increase in NH 4 observed at CB, TN concentrations were relatively stable over the same time period with no statistically significant increase in TN at this station (Fig. 5).Conversely, TN increased steadily from 2002 to 2019 at TA, with the slope of the time series trend greater than 0 throughout the study period, reflecting a gradually increasing, linear slope (Supporting Information Fig. S8).

Time series periodic components
At CB, statistically significant periodicities in NH 4 concentrations were sporadically observed on the scale of months (i.e., < 100-d periods), though based on time-averaged power spectra, CB exhibited significant signals in NH 4 only at periods of 2-3 yr at high tide and no significant periodicities at low tide (Fig. 6).Conversely, at TA, significant periodicities of 1 yr were observed for nearly every parameter at both high and low tides (Figs. 6, 7, Supporting Information Figs.S9, S10).As expected, Chl a concentration peaked on an annual scale within both estuaries during both tide heights (Supporting Information Fig. S9).At CB, DIN : DIP began exhibiting cyclic behavior at multiple periods beginning in 2010 (Supporting Information Fig. S10), but there was no evidence that periodicity systematically changed over time for any other parameter measured at CB or for any parameter at TA (Figs. 6, 7; Supporting Information Figs.S9, S10).

Drivers of ammonium across estuary types
At CB, mean annual sea level height accounted for > 90% of the importance value for mean annual NH 4 concentration at both high and low tide stages, with the multivariate ENSO index accounting for nearly all of the remainder (Fig. 8).In comparison, multiple parameters appear to be important to annual NH 4 levels at TA, and these drivers differ between high and low tides.At high tide, sea level height (66%) and MEI (32%) both contributed strongly to NH 4 based on their importance values, while at low tide, only local precipitation did not contribute at least 10% of the importance value (Fig. 8).At low tide at TA, mean flow from the Pee Dee River sharply increased in importance, sea level height decreased in importance by nearly an order of magnitude, and the multivariate ENSO index nearly doubled in importance compared with high tide conditions (Fig. 8).At CB during both tide stages, mean sea level height was strongly positively correlated with NH 4 concentrations, while at TA sea level height and NH 4 were positively correlated at high tide but did not exhibit a clear relationship at low tide (Supporting Information Fig. S11).

Discussion
Large estuaries can exhibit climate-driven signals in water quality parameters associated mainly with riverine input (Cloern et al. 1983;Boynton and Kemp 2000;Paerl et al. 2014).Here, we demonstrate that long-term increases in concentrations of some nutrient parameters and Chl a can also be observed within smaller estuarine systems dominated by oceanic inputs and that the main driver of those increases appears to be the water level within the estuary.In the minimally impacted (by local anthropogenic factors) and ocean-dominated North Inlet estuary, we document   Portions of the time series where a significant periodicity was observed are outlined (p < 0.05).Shaded regions signify the cone of influence.The bottom panels show time-averaged power spectra, with peaks representing periods that have higher power throughout the time series; the dashed red line indicates the 95% significance level based on a red-noise background spectrum.Increasing power at the longest periods is an artifact of how the width of the Morlet wavelet changes across frequencies, that is, at long periods it is narrower in frequency space and estimates higher peaks of power.Note that y-axis scales differ.
North Inlet estuary (Krask et al. 2022).Our results suggest that sea level rise can be a highly influential driver of water column NH 4 concentrations in smaller, ocean-dominated estuaries.In these systems, internal biogeochemical and hydrological dynamics between salt marshes and tidal creeks, which change in response to increasing water levels, are major controls on nutrient cycles.Conversely, NH 4 concentrations in a larger, river-dominated system (Winyah Bay) are not changing with sea level and instead appear to be largely a function of upstream watershed inputs.This is an important finding because sea level rise has previously been shown to influence salinity within a river-dominated estuary Fig. 7. Wavelet analysis to characterize periodic frequencies in total N at Clambank Landing, CB (left) and Thousand Acre, TA (right) at high (top rows for each figure type) and low tide conditions (bottom rows), from 2002 to 2019.Top panels contain continuous wavelet power spectra showing periodicities.Portions of the time series where a significant periodicity was observed are outlined (p < 0.05).Shaded regions signify the cone of influence.The bottom panels show time-averaged power spectra, with peaks representing periods that have higher power throughout the time series; the dashed red line indicates the 95% significance level based on a red-noise background spectrum.Increasing power at the longest periods is an artifact of how the width of the Morlet wavelet changes across frequencies, that is, at long periods it is narrower in frequency space and estimates higher peaks of power.Note that y-axis scales differ.
independent of changes to freshwater input rates (Ross et al. 2015).This difference is potentially due to the complex nature of the biogeochemical cycling of nutrients compared with the relatively straightforward physical properties which drive estuarine salinity gradients.Nonetheless, the interaction between sea level rise and benthic-pelagic coupling has not been the focus of substantial research, particularly within estuarine ecosystems.This is surprising considering the potential for sea level rise-related impacts as larger areas of coastal habitat experience more frequent tidal inundation.
Drivers of mean NH 4 concentration at the annual scale differed between North Inlet estuary and Winyah Bay, supporting our hypothesis that mean sea level is the dominant force leading to increases in NH 4 and Chl a concentrations within North Inlet.Previously, we demonstrated a positive relationship between the inundation of the salt marsh platform due to sea level rise and concentrations of NH 4 within the interstitial porewaters of salt marsh sediments (Krask et al. 2022).NH 4 is the dominant form of inorganic nitrogen found in salt marsh porewaters due to the highly reducing nature of the sediments, and porewater concentrations of NH 4 are often orders of magnitude higher than concentrations in adjacent tidal creeks (Morris 2000).Our measurements here, showing increases in water column NH 4 concentrations despite unchanging TN concentrations, likely reflect internal changes in biogeochemical processes that change the ratio of N constituents comprising the total pool.For example, enhanced organic matter remineralization could simultaneously drive increasing NH 4 concentrations and decreasing DON concentrations, which would result in no net change to TN concentrations.This potential explanation is supported by the differences observed between the relative proportions of monthly mean average NH 4 and DON concentrations between the 2002-2004 and 2017-2019 time periods at CB (Fig. 2).Relatedly, increasing NH 4 concentrations in the absence of increasing PO 4 concentrations, and therefore increases in the DIN : DIP ratio, could result from key differences in the biogeochemical behavior of the two nutrients.In particular, PO 4 availability is strongly influenced by the physicochemical properties of sediments, which can promote the adsorption of PO 4 (Lucotte and d'Anglejan 1988), resulting in lower concentrations.In addition, porewater flux from marsh sediments to adjacent tidal creeks is predominantly driven by tidal oscillations (reviewed by Guimond and Tamborski 2021), with models suggesting that increased tidal inundation at higher marsh elevations enhances this "tidal pumping" mechanism, leading to higher volumes draining from the marsh (Wilson and Morris 2012).Within this historically nitrogenlimited estuary (Lewitus et al. 1998;Reed et al. 2016;Pinckney et al. 2020), we suggest that increased porewater NH 4 concentrations coupled with increased porewater flux volumes, both in response to sea level rise, are leading to increases in NH 4 and Chl a in tidal creeks within North Inlet estuary.Conversely, this mechanism may be negligible with respect to NH 4 concentrations and trends over time in Winyah Bay because there is a lower ratio of intertidal marsh surface area to water volume.As such, NH 4 export equivalent to that from North Inlet per unit intertidal marsh sediment surface area is more dilute per unit volume of water in Winyah Bay.Understanding the mechanisms behind the 4-yr lag in Chl a (relative to increases in NH 4 ) is a key future research question.
Importantly, changes in other potential drivers of nutrient export to estuarine waters, such as increased precipitation, have not been observed within the North Inlet estuary.This leads us to conclude that increases in water column nutrient concentrations in this system are due to environmental changes associated with anthropogenic climate warmingnamely, sea level rise.Our analyses of potential drivers of water column nutrients at the annual scale support this hypothesis for trends in the North Inlet estuary, where sea level height contributed over 90% to the importance value driving high and low tide NH 4 concentrations.Furthermore, while we would expect that the influence of sea level rise on NH 4 is small in Winyah Bay relative to broader regional processes controlling nutrient cycling throughout the watershed, annual mean sea level height was an important contributor to NH 4 concentration at high tide when more of the marsh platform is submerged.Surprisingly, Pee Dee River flow accounted for relatively little of the importance value at TA despite being the dominant source of freshwater into Winyah Bay.We see two possible explanations for this result.First, annual mean river flow may not fully capture the dynamic nature of freshwater input into this estuary.For example, mean river flow could be identical in years characterized by regular, light precipitation (with associated low levels of continuous nutrient inputs) and more periodic but heavier rain events (with potentially different residence times or changes to hydrology).Mean annual river flow would not be reflective of these potential differences in nutrient input.Second, NO 3 is often the primary form of inorganic nitrogen entering river-dominated estuaries with watersheds strongly influenced by agriculture (Randall and Mulla 2001;Costanzo et al. 2003;Jordan et al. 2003).Thus, the minimal importance value of Pee Dee River flow on NH 4 concentrations at TA may not be reflective of the influence of river input on other nitrogen constituents (i.e., NO 3 ) or the TN pool.The increase in TN concentrations observed at TA in the absence of a concurrent increase in NH 4 indicates that the NO 3 or organic N fractions of the total pool are increasing.We document increasing NO 3 concentrations at TA (Supporting Information Fig. S7) over the course of the time series and provide evidence of consistent monthly mean average DON concentrations between the 2002-2004 and 2017-2019 time periods (Fig. 2).This likely reflects input from the large Winyah Bay watershed rather than effects of sea level rise.The effects of watershed size and composition are also evident in our analysis of periodic frequencies.Shorter frequency periods were common for many parameters at CB whereas periods of 1 yr were common at TA, suggesting conditions within the ocean-dominated North Inlet estuary are driven by shorter-term processes than in Winyah Bay, where the annual evapotranspiration cycle within the large watershed can overwhelm smaller scale, localized events.The shorter-term periodicity in NH 4 at CB could reflect an interaction between intra-annual tidal cycles and seasonal nutrient patterns.We suggest that the long-term increase in NH 4 concentrations at CB is a function of both enhanced tidal pumping (thus, greater porewater flux volumes) related to expanded tidal ranges as well as increases in marsh porewater nutrients resulting from more frequent and prolonged inundation.However, these mechanisms can occur independently of one another and may influence NH 4 concentrations to a varying degree throughout the year.Marsh porewater nutrients maintain an annual pattern closely related to seasonal changes (Morris 2000), such that average ambient porewater NH 4 concentrations are drastically lower, and the biogeochemical effect of inundation is likely minimized during cooler months.However, multiple periods of sustained increases or decreases in water level range resulting from tidal variability (i.e., perigean spring tides; Ray and Merrifield 2019) occur within a year that modulates porewater flux volumes.Thus, the degree of overlap between the occurrence of both mechanisms throughout the year may in part explain subannual periodicities in NH 4 concentration.In addition to the long-term increase in DIN : DIP in the North Inlet estuary (Supporting Information Fig. S5), the frequency of occurrence of DIN : DIP molar ratios > 25 : 1 has been steadily increasing since 2010 (unpublished data), far above the Redfield ratio of 16 : 1.Although we recognize that caution is warranted in making assertions about phytoplankton nutrient limitation based solely on inorganic N : P ratios (Ptacnik et al. 2010;Bell et al. 2018), the substantial increase in NH 4 at CB raises the possibility of a shift in degree of N vs. P limitation in this estuary (Lewitus et al. 1998;Pinckney et al. 2020).In addition, shifts in the proportion of inorganic (highly bioavailable) N and increasing DIN : DIP could have important consequences for phytoplankton biomass, production, and community composition (Glibert et al. 2016;Van Meerssche and Pinckney 2019).These changes could in turn have consequences for estuarine trophodynamics given the potential for prey selectivity by grazers on phytoplankton (Fileman et al. 2010;Siuda and Dam 2010).Finally, given that changing proportions of nutrient constituents could alter phytoplankton community structure, these changes in abundance or composition of phytoplankton could reverberate throughout food webs via a bottom-up trophic cascade (Stroud et al. 2019).The extent to which relative increases in DIN input to estuaries historically dominated by organic N shifts nutrient limitation and trophic structure is a critical yet unresolved question.
Sea level rise-induced changes to nutrient concentrations could impact the broader suite of ecosystem services and functions associated with estuaries.For example, eutrophication of relatively pristine estuaries which have minimal anthropogenic impacts within their watersheds could lead to unexpected increases in harmful algal blooms, which are a growing stressor in coastal ecosystems and can interact synergistically with climate change (Griffith and Gobler 2020).Similarly, changes to nutrient cycling could alter the growth and distribution of salt marsh plants (Morris et al. 2013) and thus impact habitat value for important fishes, invertebrates, and avifauna which rely on the structure provided by emergent vegetation within marshes.Finally, marsh drowning due to coastal squeeze (Pontee 2013) could, in theory, lead to estuarine oligotrophication based on the sequence of events we outline here (i.e., pumping of nutrients from marsh porewater into adjacent estuarine waters).However, we do not believe this to be a likely scenario given the overwhelming levels of terrestrially-derived, anthropogenic nutrients simultaneously entering estuaries.In either case, the ability to identify these ecological impacts is dependent on long-term time series data (Lindenmayer et al. 2012), and this study demonstrates the value of rigorous monitoring programs for detecting and describing environmental changes associated with sea level rise as well as anticipating how coastal ecosystems will respond to changing climatic and hydrologic regimes.

Fig. 1 .
Fig. 1.Map of focal estuaries: North Inlet and Winyah Bay, located in Georgetown County, South Carolina, USA.The first panel shows the entire Winyah Bay watershed; the red box highlights the focal estuaries.The second panel shows Winyah Bay proper and the entire North Inlet estuary watershed, as well as the monitoring stations at Clambank Landing and Thousand Acre.Sea level height and precipitation data were collected at Oyster Landing, also denoted.Imagery credits: Esri, USDA Farm Service Agency, HERE, Garmin, Food and Agriculture Organization, NOAA, USGS, EPA, and National Park Service.This map was created using ArcGIS ® software by Esri.ArcGIS ® and ArcMap™ are the intellectual property of Esri and are used herein under license.Copyright © Esri.All rights reserved.

Fig. 2 .
Fig. 2. Comparison of quantities of TN constituents between North Inlet estuary (Clambank Landing station, CB) and Winyah Bay (Thousand Acre station, TA) during the initial 3 yr (2002-2004) and final 3 yr (2017-2019) of the study.Bars show monthly mean values.TN, TDN (total dissolved N-both inorganic and organic), NO x , and NH 4 (dissolved inorganic N) were all measured directly.TDN is measured by the same analysis described in the text for TN, only on a filtered sample.Thus, DON = TDN À (NH 4 + NO x ) while PN = TN À TDN.

Fig. 3 .
Fig. 3. Partial smooths, time series, and change-point analysis for GAM of NH 4 at Clambank Landing low tide sampling events.(A, B) Partial effect plots show s(day of time series) and s(month) with shading representing confidence intervals (CI) and associated residuals.Numbers on y-axes labels of partial effect plots represent the estimated degrees of freedom for that particular smooth term.In panel (A), "Day" corresponds to the number of days since 01 January 1970, with the first data point in the time series on day 11,693 (equivalent to 06 January 2002).(C) Fitted values from GAM are overlaid on the raw time series data.(D) Change-point analysis shows the first derivative and simultaneous CIs on the estimated trend for the time series (after removing the effect of the monthly smooth term).Locations where the CIs do not include 0 (horizontal black line) indicate periods of significant change in the trend.

Fig. 4 .
Fig. 4. Partial smooths, time series, and change-point analysis for GAM of Chl a at Clambank Landing low tide sampling events.(A, B) Partial effect plots show s(day of time series) and s(month) with shading representing confidence intervals (CI) and associated residuals.Numbers on y-axes labels of partial effect plots represent the estimated degrees of freedom for that particular smooth term.In panel (A), "Day" corresponds to the number of days since 01 January 1970, with the first data point in the time series on day 11,693 (equivalent to 06 January 2002).(C) Fitted values from GAM are overlaid on the raw time series data.(D) Change-point analysis shows the first derivative and simultaneous CIs on the estimated trend for the time series (after removing the effect of the monthly smooth term).Locations where the CIs do not include 0 (horizontal black line) indicate periods of significant change in the trend.

Fig. 5 .
Fig. 5. Partial smooths, time series, and change-point analysis for GAM of total N at Clambank Landing low tide sampling events.(A, B) Partial effect plots show s(day of time series) and s(month) with shading representing confidence intervals (CI) and associated residuals.Numbers on y-axes labels of partial effect plots represent the estimated degrees of freedom for that particular smooth term.In panel (A), "Day" corresponds to the number of days since 01 January 1970, with the first data point in the time series on day 11,693 (equivalent to 06 January 2002).(C) Fitted values from GAM are overlaid on the raw time series data.(D) Change-point analysis shows the first derivative and simultaneous CIs on the estimated trend for the time series (after removing the effect of the monthly smooth term).Locations where the CIs do not include 0 (horizontal black line) indicate periods of significant change in the trend.
increases in NH 4 and Chl a from 2002 to 2019, although we observed no change in TN concentration.In comparison, in the river-dominated Winyah Bay, NH 4 and Chl a concentrations exhibited higher variability throughout the study period but have not increased over time.However, TN concentrations have increased steadily in Winyah Bay, indicating increased input of a nitrogen constituent other than NH 4 to the system, and in fact, we document increases in NO 3 concentrations at our Winyah Bay sampling station but not within North Inlet.Sea level rise is driving higher water levels in estuaries across the southeast U.S.(Sweet et al. 2022), including the

Fig. 6 .
Fig. 6.Wavelet analysis to characterize periodic frequencies in NH 4 at Clambank Landing, CB (left) and Thousand Acre, TA (right) at high (top rows for each figure type) and low tide conditions (bottom rows), from 2002 to 2019.Top panels contain continuous wavelet power spectra showing periodicities.Portions of the time series where a significant periodicity was observed are outlined (p < 0.05).Shaded regions signify the cone of influence.The bottom panels show time-averaged power spectra, with peaks representing periods that have higher power throughout the time series; the dashed red line indicates the 95% significance level based on a red-noise background spectrum.Increasing power at the longest periods is an artifact of how the width of the Morlet wavelet changes across frequencies, that is, at long periods it is narrower in frequency space and estimates higher peaks of power.Note that y-axis scales differ.

Fig. 8 .
Fig. 8. Random forest output comparing importance values (mean + SE) for potential drivers of NH 4 at Clambank Landing (A, B) and Thousand Acre (C, D) at high tide (gray fill) and low tide (white fill), respectively.Mean importance values and associated standard errors were calculated from 1000 runs for each site/tide level combination.

Table 1 .
Generalized additive model output for models fit to NH 4 , NO 3 , DIN : DIP, Chl a, and TN time series from Clambank Landing and Thousand Acre, within the North Inlet-Winyah Bay National Estuarine Research Reserve.Models fit to data collected at high tide are shown in gray and low tide models are in white.Estimated degrees of freedom (edf) is a measure of the "wiggliness" of the partial effect of a given model term.