Elevated organic carbon pulses persist in estuarine environment after major storm events

Abstract

Estuaries regulate transport of dissolved organic carbon (DOC) from land to ocean. Export of terrestrial DOC from coastal watersheds is exacerbated by increasing major rainfall and storm events and human activities, leading to pulses of DOC that are shunted through rivers downstream to estuaries. Despite an upward trend of extreme events, the fate of the pulsed terrestrial DOC in estuaries remains unclear. We analyzed the effects of seven major tropical cyclones (TC) from 1999 to 2017 on the quantity and fate of DOC in the Neuse River Estuary (NC, USA). Significant TC‐induced increases in DOC were observed throughout the estuary; the increase lasting from around 50 d at head‐of‐tide to over 6 months in lower estuary. Our results suggest that pulsed terrestrial DOC associated with TCs temporarily overwhelms the estuarine filter's abiotic and biotic degradation capacity under such high flow events, enhancing the shunt of terrestrial carbon to the coastal ocean.

---

Scientific Significance Statement

The recently formulated “new normal” with more frequent and more intensive precipitation associated with tropical cyclones, leads to high loads of terrestrial dissolved organic carbon (DOC) into coastal ecosystems. A part of this DOC is very bioavailable, causing higher pCO2 in the estuarine waters that persists for weeks following a storm event. However, there are major knowledge gaps about the fate of the pulsed DOC, one of the major obstacles being the lack of locations with sufficient water quality monitoring data available. Our key result is that the elevated DOC lasts for 50 ‐200 days in the estuary, and is controlled mainly by mixing with seawater − not by biogeochemical processing − thus overwhelming abiotic and biotic degradation capacity of the estuarine filter under high flow events.

---

Estuaries are important ecosystems that mediate processing of terrestrial carbon (C) at the land‐ocean interface (Massicotte et al. 2017). CO₂ emissions associated with estuarine organic carbon (OC) mineralization and sedimentary burial are major components of the global C cycle (Chen et al. 2012; Bianchi et al. 2018). OC inputs from catchments to estuaries are largely controlled by hydrology (precipitation and consequent river discharge) and human watershed activities; however, major storms and tropical cyclones (TCs) can cause marked increases in OC loading to estuaries compared to nonstorm years (Hinton et al. 1997). Especially “wet” TCs that provide greater than the 90^th percentile of weekly average river flows for the Atlantic storm season (May–October) cause considerable increases in OC concentrations (Paerl et al. 2018). These low‐frequency, high‐impact events are responsible for a significant percentage of annual terrestrial DOC input to drainage networks (Raymond et al. 2016).

Climate change is altering hydrological cycles globally; for instance, in coastal Mid‐Atlantic US a period of unprecedented TC‐related high precipitation has been prevailing since the late 1990s (Paerl et al. 2019). The high overall precipitation is linked with a trend toward higher frequency of TCs and intense precipitation associated with them (Webster et al. 2005; Wuebbles et al. 2014), which represents a recent regime shift in coastal hydrology with major ramifications for carbon cycling. TCs change hydrological flow paths in the catchment, causing significant and rapid increases in river network stream velocities, and shunting large quantities of pulsed DOC downstream (Buffam et al. 2001; Raymond et al. 2016). The pulse‐shunt episodes increase the quantity and alter the bulk composition of DOC drained to estuaries, introducing DOC that is fresher, i.e., less diagenetically altered (Sobczak and Findlay 2002; Peter et al. 2012; Hitchcock and Mitrovic 2015; Rudolph et al. 2020). In particular, the flooding of coastal wetlands during storm events causes drastic changes in carbon biogeochemistry as large pools of labile DOC are released from wetlands to estuaries (Osburn et al. 2019a).

Under baseflow conditions, up to 20% of the changes in estuarine DOC characteristics is estimated to occur through biogeochemical processing within the estuary (Asmala et al. 2016). During pulse events, the effects of storm‐driven DOC inputs extend from estuaries to adjacent coastal waters, where pulsed DOC causes higher mineralization rates lasting weeks after the storm (Osburn et al. 2019b). Overall, the fate of the pulsed terrestrial DOC in the receiving estuarine and coastal environments is unclear. Due to high internal spatial and temporal heterogeneity in the biogeochemical processing of DOC, it is difficult to quantify the carbon balance of even a single estuary (Bauer et al. 2013). A key challenge is to delineate the elevated DOC concentrations from background concentrations. However, long‐term time series observations across estuarine salinity gradients can be used to quantify deviations from the expected baseline values (Asmala et al. 2018).

We used 2410 observations during 1999–2017 on the Neuse River Estuary (NRE) (North Carolina, USA) to quantify the effect of six major TCs on elevation of its DOC concentrations. Furthermore, we used a breakpoint analysis to resolve the time needed for the system to return to the baseline DOC levels. Our aim was to clarify the space–time scales over which to consider TC impacts on estuarine and downstream coastal ecosystems.

Material and methods

Study area

We obtained data analyzed here from a long‐term monitoring program, ModMon, in Neuse River Estuary, North Carolina, USA (NRE; Fig. 1). ModMon is a collaborative monitoring program between the University of North Carolina Chapel Hill and the State of North Carolina's Department of Environmental Quality (http://paerllab.web.unc.edu/projects/modmon/). Sampling was carried out by bimonthly visits to six mid‐river stations along the estuarine portion of NRE (Fig. 1), including vertical profiles and collection of near‐surface water for water quality analysis (Paerl et al. 2018).

1 Fig.. Map of the sampling sites in Neuse River Estuary (NRE0–NRE160). The catchment of Neuse River is shaded light gray.

Field and laboratory analyses

Salinity was measured using a YSI 6600 multiparameter water quality sonde (Yellow Springs). Sondes were calibrated prior to each sampling trip according to the manufacturer guidelines. Water samples were collected at 0.2 m depth using a nondestructive diaphragm pump, dispensed into 4 liter polyethylene bottles. Samples were returned to the laboratory for processing within 4 h of collection. Samples for dissolved organic carbon (DOC) analysis were filtered through precombusted (4 h at 525°C) 25 mm diameter Whatman GF/F filters into acid‐washed and sample‐rinsed 150 mL polyethylene bottles that were subsequently frozen at −20°C. DOC was measured on a Shimadzu Total Organic Carbon Analyzer (TOC‐5000A). The dataset consisted of 2410 DOC observations in surface waters. A summary of the salinity and DOC observations is presented in Table S1. Neuse River discharge from the USGS gauge site at Ft. Barnwell (Station ID: 02091814) was available from October 2000 onwards (thus not covering the study period), and was used to calculate DOC fluxes by multiplying daily average discharge with the measured DOC concentration (Fig. S1). Due to incomplete discharge data and additional uncertainties introduced with loading calculations, DOC loading data was not included in the main analysis.

Statistical analyses and baseline calculation

We tested the DOC time series for a potential monotonic trend with a seasonal Kendall test, which accounts for seasonality by computing the Mann–Kendall test on each month separately. A nonparametric modeling approach, generalized additive model (GAM), was used to explore the possible seasonality in DOC observations. GAM was carried out using the gam function in the mgcv package for R (Wood and Wood 2015), which solves the smoothing parameter estimation by using the generalized cross validation (GCV) criterion. Day of the year was used as the explanatory variable and DOC concentration as the response variable, the number of smoothing terms (splines) was chosen by GCV in the GAM procedure, and cyclic cubic regression splines were used for smoothing across the annual cycle. After presence of seasonality was confirmed with the GAM approach, years with and without TCs were separated and analyzed for differences in the annual cycle (Fig. 2). An ANOVA was used to test differences among months between years with and without TCs. A baseline for DOC concentrations in NRE was calculated for each station and month by using mean values from 1999 to 2017 (Fig. S2), excluding years with wet TCs (final number of DOC observations = 14–27 per group). The difference between the monthly mean and median DOC concentrations among months for each station ranged from −20 to 90 μmol L⁻¹ (equaling −3.8 to 16.7% of total DOC concentration), and the average difference was 22 μmol L⁻¹ (3.9%). Further, Shapiro–Wilk test of normality shows that out of 72 Station‐Month observation groups, 42 were normally distributed (Fig. S3). The relatively small difference between mean and median values indicates that extremely high (or low) observations were not common, and do not disproportionally affect the baseline calculations. It should be noted that the determination of a DOC baseline for highly dynamic systems such as estuaries inevitably contains large uncertainties resulting from, e.g., unaccounted storm runoff related diffuse sources (Rudolph et al. 2020) or heterogeneity in riverborne water parcels (Massicotte and Frenette 2013). The calculated values were used to quantify the DOC deviations from baseline concentrations after TCs. We used segmented regression analysis to identify the breakpoints, i.e., the intercept of the piecewise linear DOC deviation and the long‐term baseline (Fig. 2). In this approach, the value of DOC deviation (y) was assumed zero when time (x; number of days since last storm) became independent, and the slope of the linear regression was set to be zero on the right‐hand side of the breakpoint (the unknown values). A segmented regression model was carried out with the segmented package for R (Muggeo 2008). Simple linear regression was used for analysis of relationships between the return to baseline, salinity and distance from head‐of‐tide. All statistical analyses were conducted using R Software (R Core Team 2020).

2 Fig.. Interannual and seasonal variations of dissolved organic carbon (DOC) concentration in freshwater end‐member of Neuse River Estuary (NRE 0). (a) Time series from 1999 to 2017, with wet TCs (sensu Paerl et al. 2018) is indicated with vertical lines. (b) Seasonal dynamics in the NRE 0 contrasting years with wet TCs (solid blue circles) and years without wet TCs (white circles). Presence of seasonal trends was confirmed with a generalized additive model, which showed that the seasonal component was significant (p < 0.001).

Results

Temporal DOC dynamics

Overall, there was a significant monotonic positive trend in the DOC concentrations at all station except in the most downstream NRE 160 (an annual increase of 3.8–6.8 μmol L⁻¹ in the mean concentration; Fig. S4). We observed large variability in DOC concentrations

μmol L⁻¹, ranging between 214 and 1543 μmol L⁻¹ (Fig. 2a). During the 1999–2017 study period, seven major wet TCs occurred in the region (Fig. 2; Paerl et al. 2018). All occurred within a 6‐week time window during the annual cycle between August 27 and October 8. Following each TC, except for a more coastal storm, Gordon in 2000, there was an apparent increase in Neuse River DOC concentration (Fig. 2), as well as in DOC loading (Fig. S1). Detailed analysis of DOC loadings instead of concentrations might yield additional insight about estuarine carbon cycling, but due to data limitations it was not possible in this study. The observed DOC increase was transient, and DOC concentrations returned to prestorm baseline levels shortly after the event. The annual cycle for DOC showed different patterns for years with and without TCs (Fig. 2b). In years without TCs, the highest DOC concentrations were observed in late summer (August–September), with the lowest concentrations in winter (December–January). The elevated DOC concentrations due to TCs were most apparent in autumn (September–November), when the difference between nonstorm years and storm years was ~ 200 μmol C L⁻¹. Post hoc comparisons using the Tukey HSD test indicated that the monthly mean DOC concentration for years with TCs was significantly (p < 0.02) different from years with no TC during the most active September and October months (Fig. S5). Differences were nonsignificant (p > 0.05) for other months. Mean (± 1 SD) DOC concentrations during wet TC years were 910 ± 282 and 809 ± 254 μmol L⁻¹ for September and October, respectively. In years without wet TCs, DOC concentrations were 669 ± 192 and 577 ± 137 μmol L⁻¹ in September and October, respectively.

DOC return to baseline

To quantify the increase in DOC concentrations in the estuary after each TC, we pooled the observations following seven wet TCs (sensu Paerl et al. 2018) and quantified the deviations of DOC concentration from the long‐term baseline. Segmented regression analysis revealed that these deviations follow a pattern where initially high deviations gradually decrease linearly and reach the baseline (“breakpoint”) DOC concentration in the estuary (Fig. 3). The intercept value of the segmented regression provides an estimate of the mean deviation of DOC from the baseline at t₀, i.e., at the time when the TC made landfall in the catchment and caused the DOC pulse in the NRE. The initial DOC deviations were highest at head‐of‐tide (NRE 0; Fig. 3a), on average 532 μmol C L⁻¹ higher than the baseline concentration. The deviations at t₀ decreased gradually down‐estuary, and were 270 μmol C L⁻¹ in the lower estuary, where it flows into Pamlico Sound (Fig. 3f). We used the breakpoint of the segmented regression analysis to determine the time required for the DOC concentrations to return to baseline concentrations after TCs. After reaching the breakpoint, DOC concentration deviations vary around zero, indicating a return to the long‐term baseline level. At head‐of‐tide (NRE 0), the effect of the TC‐induced DOC pulse diminished on average after 55 d (Fig. 3a). The breakpoint time to return to baseline DOC concentrations increased down‐estuary, being 205 d near Pamlico Sound (NRE 160; Fig. 3f). Thus, the decreasing intercept and slope values of the regression lines in a seaward direction show a smaller initial deviation and longer lingering of the pulsed DOC in the lower estuary.

3 Fig.. Segmented regression analysis of estuarine DOC deviations after seven major TCs between 1999 and 2017 along Neuse River Estuary (a–f). Points indicate deviations from the baseline concentration, which is established separately for each sampling station and month. Segmented regression analysis (red solid line) is used to identify the breakpoints, which indicate the return of the DOC concentrations to the baseline after storm events. Breakpoint value (in days) and standard error (gray shaded area) are given for each sampling station, along with coefficient of determination and formula for the first segment of each segmented linear regression.

Relationship between lingering DOC and hydrology

We used the information from the segmented regression analysis to investigate potential underlying mechanisms, and found two strong linear relationships between the time it takes the DOC concentrations to return to the baseline levels and (1) the station long‐term mean salinity (Fig. 4a), as well as (2) the station's distance from head‐of‐tide where DOC loads originate (Fig. 4b). Here, mean salinity and distance from the head‐of‐tide can be considered as proxies for hydrological distance between various locations in the system. Conceptually, the linear relationship between the time to baseline return and hydrological distance suggests a ratio between the rate of DOC degradation (D) and the rate of the DOC supply (S) (D : S) is very close to zero (Fig. 4c). If biogeochemical processing (e.g., biodegradation, photodegradation) significantly decreased the amount of pulsed DOC during its passage through the estuary (resulting in D : S > 0), the resulting response curve would be nonlinear. In general, the slope value of D : S is likely unique for individual estuaries based on local hydrogeomorphology, but in general a linear relationship is indicative of a pulse‐shunt type process (Hounshell et al. 2019; Rudolph et al. 2020). As an estuary deepens and widens toward the coast, the initially high concentrations due to pulsed riverine DOC are diluted and dissipated due to the hydrological forcing as the pulse is attenuated down‐estuary (Fig. 4c inset).

4 Fig.. The relationship between baseline return period length and estuarine hydrology. (a) Return to baseline and long‐term mean salinity of each sampling station. (b) Return to baseline and distance from freshwater end‐member of the NRE. (c) Conceptual schematic, depicting anticipated relationships between baseline return and hydrological distance, given different ratios between the rate of DOC degradation and the rate of DOC supply (D : S). Inset figure illustrates a time sequence of the passage of a DOC (yellow) pulse through an estuary, during which the initially high concentrations are diluted and dissipated during the transit along the estuarine gradient.

Discussion

We observed an upward decadal trend in the DOC concentrations in the Neuse River Estuary, which concurs with recent findings of increasing organic matter concentrations throughout the northern hemisphere (Filella and Rodríguez‐Murillo 2014). Drivers for increasing DOC concentrations may be both long‐term (e.g., climate change driven DOC release from soils; Hongve et al. 2004; Evans et al. 2006; Asmala et al. 2019) and short‐term (e.g., weather driven alteration of hydrological flowpaths in heterogeneous catchments; Laudon et al. 2011). The main driver of these increases in DOC concentration likely is hydrological connectivity of the river channel to adjacent wetlands, which hold vast quantities of organic matter (Majidzadeh et al. 2017; Rudolph et al. 2020). Exacerbating the general increase in DOC concentrations, storm events are being recognized as important drivers of organic carbon loading from catchments to drainage networks, but the magnitude of DOC delivery to coastal environments attributable to these events remains uncertain (Avery Jr et al. 2004; Clark et al. 2007; Brown et al. 2014; Paerl et al. 2018). Our results suggest that elevated DOC concentrations after major storms are attenuated mostly due to dilution of river water by mixing with sea water, and not by biological degradation or physicochemical removal processes.

While the estuarine biological filter performs substantial organic matter processing and mineralization with respect to quantity and complexity (Smith and Hollibaugh 1993; Montgomery et al. 2013; Osburn et al. 2019a), it is unable to significantly process the surge of DOC associated with TCs, overwhelming the estuarine biogeochemical filter. Compared to normal baseline flow conditions, the increased flow associated to TCs result in decreased residence time within the estuary, explaining the exceedance of the estuarine filter. The linear relationship between hydrological distance and time to baseline return indicates that the net removal rate is far too low to cause a marked decrease in the DOC pool, as rate of supply vastly increases rate of degradation (D : S ≈ 0). This outcome overwhelms the estuarine filter and leads to downstream processing or storage in the coastal or open ocean, rather than estuarine DOC degradation and CO₂ generation, at least until freshwater discharge returns to baseline conditions. Under those conditions, residence time in the estuary is long enough to allow microorganisms to exploit increased DOC concentrations and photodegradation to occur. Hence, the surge of water and DOC to estuaries outpaces the biogeochemical processing of DOC in these systems well after the storm has passed (Osburn et al. 2019b).

Our observations and segmented modeling approach provide a mechanistic basis for estuarine biogeochemical observations related to storm events, including elevated DOC concentrations (Paerl et al. 2018), altered DOC characteristics (Fellman et al. 2009; Nguyen et al. 2010) and higher CO₂ degassing (Mead and Wiegner 2010; Crosswell et al. 2014; Osburn et al. 2019b). We modify the view that major TC events profoundly alter the carbon biogeochemistry of the estuaries temporarily, as DOC from the high‐discharge events enters a passive pipe and is transported conservatively and shunted through the estuary. Thus, the delivery of DOC to the coastal ocean during these transient high flow events ultimately is set by loadings at the head of tide, and remains under control by simple physical mixing in response to changing residence time. Because this study focused only on the dissolved fraction, the transport and fate of particulate organic carbon associated to TCs remains poorly constrained. Further, accurate mass flux calculations would allow a strong quantitative estimate of how much DOC actually is delivered to the estuary, degraded within the estuary, and/or exported from the estuary to the coastal ocean. However, monitoring the response of estuaries as outlined in this work provides the basis for understanding how increasing intensity of TCs can alter coastal carbon cycling.

Acknowledgments

We appreciate the assistance of J. Braddy, B. Abare, R. Sloup, students, and technicians that participated in field and laboratory work supporting this publication. This research was funded by Academy of Finland Project 309748, US National Science Foundation Projects 0825466, 0812913, 1338606, 1705972, 1706009, The North Carolina Dept. of Environmental Quality and National Fish and Wildlife Federation (Project 53916: ModMon Program), the Lower Neuse Basin Association, NC Sea Grant, and the UNC Water Resources Research Institute.