Use of Geospatial, Hydrologic, and Geochemical Modeling to Determine the Influence of Wetland-Derived Organic Matter in Coastal Waters in Response to Extreme Weather Events

Flooding from extreme weather events (EWE), such as hurricanes, exports large amounts of dissolved organic matter (DOM) to both estuaries and coastal waters globally. Hydrologic connectivity of wetlands to adjacent river channels during flood events can potentially have a major control on the DOM exported to coastal waters after EWEs. In this study, a geographic information system based flood model was used to: (1) determine the volume of flooded wetlands in a river corridor following Hurricane Matthew in 2016; (2) compute the resulting volume fluxes of DOM to the Neuse River Estuary-Pamlico Sound (NRE-PS), in eastern North Carolina and (3) use the flood model to quantify the wetland contribution to DOM export. The flood model-derived contributions were validated with a Bayesian Monte Carlo mixing model combining measurements of DOM quality: specific UV Absorbance at 254 nm (SUVA254), spectral slope ratio (SR), and stable isotope ratios of dissolved organic carbon (δ13C-DOC). Results indicated that (1) hydrologic connectivity of the freshwater riparian wetlands caused the wetlands to become the primary source of organic matter (OM) that was exported into the NRE-PS after Matthew and (2) this source lingered in these coastal waters in the months after the storm. Thus, in consideration of the pulse-shunt concept, EWE such as Hurricane Matthew cause pulses of DOM from wetlands, which were the primary source of the OM shunted from the terrestrial environment to the estuary and sound. Wetlands constituted ca. 48% of the annual loading of DOC into the NRE and 16% of DOC loading into the PS over a period of 30 days after Hurricane Matthew. Results were consistent with prior studies in this system, and other coastal ecosystems, that attributed a high reactivity of DOM as the underlying reason for large CO2 releases following EWE. Adapting the pulse-shunt concept to estuaries requires the addition of a “processing” step to account for the DOM to CO2 dynamics, thus a new pulse-shunt process is proposed to incorporate coastal waters. Our results suggest that with increasing frequency and intensity of EWE, strengthening of the lateral transfer of DOM from land to ocean will occur and has the potential to greatly impact coastal carbon cycling.


INTRODUCTION
The frequency of extreme weather events (EWE), such as tropical cyclones, has increased in the southeastern United States (US) in recent decades (Paerl et al., 2018) and since the mid-1990s the coast of North Carolina has experienced many such events. Paerl et al. (2018) classified EWEs occurring over the past 20 years in eastern NC into three categories (dry and windy, wet and calm, and wet and windy) and compared each category to non-storm conditions to show how different types of storms would affect the Neuse River Estuary (NRE), the main tributary to Pamlico Sound (PS), 2nd largest estuarine complex in the United States. Hurricane Matthew, which occurred in October 2016, struck eastern North Carolina and its coastal waters as a Category 1 storm causing ca. $1.6 billion in damage and 28 deaths in North Carolina alone (Stewart, 2017). In their analysis, Paerl et al. (2018) found that "wet and windy" storms, such as Matthew, showed statistical differences in carbon (C) and inorganic nitrogen concentrations, partial pressure of carbon dioxide (pCO 2 ), and pH in the NRE, indicating such EWEs have major impacts on estuarine and coastal water biogeochemistry in the weeks to months after the storms.
Precipitation from EWE is projected to rise by ca. 7% with each degree Celsius of warming in the future (Prein et al., 2017). Flooding from heavy precipitation associated with EWE can result in a flushing of terrestrial dissolved organic matter (DOM) from land to coastal waters where substantial biogeochemical processing can occur returning CO 2 to the atmosphere . Depending on the residence time of floodwater in receiving waters such as estuaries and large river plumes, DOM degradation can persist for weeks to months following these events (Osburn et al., 2019a). Thus, these episodic large storm events can be significant, yet poorly constrained, influences on coastal C budgets, via terrestrial DOM's transport to coastal waters (Raymond et al., 2016).
Abbreviations: δ 13 C-DOC, stable isotope ratio of dissolved organic carbon ( ); [DOC], concentration of dissolved organic carbon (mg L −1 ); A T , total flooded area; A W , area of flooded wetlands; BMC, bayesian monte carlo; C FB , DOC concentration from river (mg L −1 ); C W , DOC concentration from wetland (mg L −1 ); DEM, digital elevation model; EWE, extreme weather events; f o , fractional contribution of marine sourced organic matter to BMC mixing model; f r , fractional contribution of river sourced organic matter to BMC mixing model; f w , fractional contribution of wetland sourced organic matter to BMC mixing model; L DOCR , dissolved organic carbon mass loading from riverine sources (Gg C); L DOCT , total dissolved organic carbon mass loading from wetland and riverine sources combined (Gg C); L DOCW , dissolved organic carbon mass loading from wetland sources (Gg C); NR, neuse river; NRE, neuse river estuary; PS, pamlico sound; PSC, pulse-shunt concept; REM, relative elevation model; S R , spectral slope ratio; SUVA 254 , specific UV absorbance at 254 nm (L mg C −1 m −1 ); V nW , volume of flooded non-wetland area; V T , total flooded volume; V W , volume of flooded wetlands.
For DOM, changes to its sources following EWEs can influence its biogeochemical processing in coastal waters (Osburn et al., 2012;Bianchi et al., 2013). In part, DOM's importance in the coastal C cycle arises from its central role in providing a basal resource to the "microbial loop". Results from Hurricane Floyd in 1999 indicated increased microbial respiration following floodwaters from the Cape Fear River in eastern NC exported into Onslow Bay (Avery et al., 2004). Bianchi et al. (2013) found a consistent post-flood result in shelf waters of the northern Gulf of Mexico, and suggested that soil-sourced DOM was rapidly lost to the atmosphere as CO 2 via photochemical and microbial processing. As a result of high rates of carbon sequestration and sediment trapping, wetlands are major C sinks that store between 350-535 Gt C (Gorham, 1995;IPCC, 2001;Mitra et al., 2005). Wetlands likely play major role in the passive pipe scenario used to describe rivers in the pulseshunt concept (PSC; Raymond et al., 2016), which builds on earlier concepts of how floods can influence biogeochemistry (Vannote et al., 1980;Junk et al., 1989). While the pulseshunt concept was initially only applied to upland streams and did not account for wetlands, estuaries or coastal waters, recent work suggests that estuaries may act in similar fashion . Delineating between wetland and upland DOM will be important to evaluate the reactivity and cycling of DOM transported to coastal waters following EWEs. Here we posit that flooding from EWEs mobilizes DOM stored in wetlands into coastal waters where it can be processed back into CO 2 .
The complexity of doing so requires a means to assess when hydrologic-connectivity occurs between water from a river and water from wetlands, via either surface waters, groundwater, or both (Pringle, 2003). During floods, an isolated wetland's hydrological-connection to a river's main channel allows the flushing of its accumulated recalcitrant OM and nutrients (Tockner et al., 1999;Wolf et al., 2013). However, recalcitrant OM in one system may become labile in another system (Boyd and Osburn, 2004); microorganisms in coastal waters are poised to break down the recalcitrant wetland OM to CO 2 , returning stores of terrestrial carbon to the atmosphere (Bauer et al., 2013).
Because wetlands are important reservoirs of OM in coastal watersheds that can become OM sources following EWEs, their contribution to the terrestrial DOM flux must therefore be determined. Osburn et al. (2019a) used geochemical proxies for organic matter (light absorption and stable C isotopes) to posit that wetlands constituted a significant fraction of the terrestrial DOM exported into Pamlico Sound via the Neuse River Estuary following Hurricane Matthew. The main objective of this study was to validate their supposition with a quantitative model. We hypothesized that the hydrologic-connectivity of the lower Neuse River to its adjacent riparian wetlands resulted in wetlands as a dominant source of DOM transported in floodwaters from extreme precipitation events, such as Hurricane Matthew. Here, the hydrologic-connectivity of riparian wetlands to a major coastal river in eastern NC was determined by a relative elevation flood model used to determine the extent of flooding that occurred in response to Hurricane Matthew in 2016. Fluxes of DOM to the NRE from its watershed were determined using concentrations of DOC from the river and wetlands and volumetric estimates of each of these two sources. This geospatial approach to compute DOM sources was validated by comparison to a Bayesian Monte Carlo mixing model (Arendt et al., 2015). Sources of OM in this system were determined with bulk carbon stable isotopes of DOC (δ 13 C-DOC), specific UV absorbance at 254 nm (SUVA 254 ), and spectral slope ratio (S R ), analyses that characterize different fractions of OM and have varying rates of source specificity (Bianchi and Canuel, 2011). By implementing the coupled geospatial-geochemical model and these analyses, the change in quality and quantity of OM from the wetlands to the NRE-PS resulting from Matthew's floodwaters was quantified. Previous studies of biogeochemical cycling and past extreme events in the NRE-PS provided important context for this study (e.g., Paerl et al., 2018;Osburn et al., 2019a).

Study Site and Sample Collection
This study was conducted in the main-stem of the Neuse River in the Coastal Plain of North Carolina, United States (Figure 1; Osburn et al., 2019a). Along the lower reach of the river above the head of tides, abundant freshwater riparian wetlands fringe the main channel of the river, before it flows southeast into the NRE-PS (NC-DEQ, 2018). The Neuse River Estuary is the largest tributary to the Pamlico Sound (Paerl et al., 2018). The Pamlico Sound is the larger portion of the Albemarle-Pamlico Estuarine System (APES), which is the second largest estuarine complex in the continental United States (NC-DEQ, 2018). Characteristics of each system are summarized in Supplementary Table S1. Both the NRE and PS are relatively shallow (average depths are ca. 3 and 5 m, respectively). There is little tidal effect (0.3-0.6 m; Giese et al., 1985) in this microtidal estuarine complex due to embayment by the barrier islands off the coast of NC (NC-DEQ, 2018). Mixing in the system is strongly affected by winds (Dixon et al., 2014). Moreover, the Neuse River's discharge can vary greatly over the course of a year (ca. 17 to 1400 m 3 s −1 for 2016). The embayment of this lagoonal system under low to modest river flow causes the NRE-PS to generally have long water residence times; however, during high flow events like Hurricane Matthew, the residence time can decrease greatly (Supplementary Table S1).
Access to a long-term geochemical dataset from routine monitoring of the NRE-PS (the Neuse River Estuary Modeling and Monitoring Project, "ModMon") provided a long-term record against which we compare the system's response to the storm. Surface water samples were collected at sites along the main axis of the Neuse River and NRE-PS over a period of 3 months, between October and December 2016, following the passage of Hurricane Matthew on October 9, 2016 (Figures 1A,B). Sampling was conducted ca. weekly across the NR, NRE, and PS sites (Supplementary Table S2). Using a bucket riverine samples were collected in 1 L brown HDPE bottles, from bridge overpasses along the main-stem of the NR at United States Geological Survey (USGS) gauged locations (Clayton, Smithfield, Goldsboro, Kinston, and Fort Barnwell; Figure 1B).
Estuarine and sound samples were collected during ModMon surveys of the NRE-PS by pumping water using a peristaltic pump from the surface, stored in 1 L opaque HDPE bottles and frozen for shipment to North Carolina State University ( Figure 1B). Collection bottles were cleaned in a detergent bath, rinsed profusely with Milli-Q ultrapure water, and left to air dry before sampling. Additional surface water sampling was completed in an area of freshwater riparian wetlands between the last riverine site on the NR near Fort Barnwell, NC, United States and head-of-tide for the estuary at the Street's Ferry Bridge crossing near Vanceboro, NC, United States in March and October 2017 ( Figure 1C).

Sample Processing and Optical Analyses
Surface water samples were thawed at room temperature, and a measured volume of water was filtered through pre-combusted (at 450 • C for 5 h) 0.7 µm Whatman glass fiber filters (GF/F) via vacuum. Prior to sample filtration, 150 mL of Milli-Q water was used to rinse each filter. The filtrate was collected into 60 mL polycarbonate bottles (detergent-washed and rinsed thoroughly with Milli-Q water) for optical analyses and 40 mL amber-tinted borosilicate glass vials (detergent-washed, rinsed thoroughly with Milli-Q water, and combusted at 450 • C for 5 h) for DOC concentration [(DOC)] and stable carbon isotope (δ 13 C-DOC, see below) analyses. Filtrate was stored at 4 • C until optical measurements were made, generally within 48 h of thaw. Filtrate for DOC analysis was acidified to a pH of 2 with 85% phosphoric acid (H 3 PO 4 ), then stored at 4 • C until measurements were completed.
Absorbance measurements were made with a Varian Cary 300UV spectrophotometer in 1 cm quartz cuvettes, over a range of wavelengths (200-800 nm), and then blank corrected using Milli-Q ultrapure water. The Napierian absorption coefficient (a) was determined using the formula: where a is absorption coefficient (m −1 ), λ is the wavelength (nm), A is the absorbance measured on the spectrophotometer, and L is pathlength of the quartz cuvette (m) (Kirk, 1994;Osburn and Morris, 2003). Helms et al. (2008) suggested that an absorption coefficient at 254 nm (a 254 ) could be used as a proxy for DOC concentration, due to high correlation of DOC concentration and a 254 . The spectral slope ratio (S R ) was calculated by dividing the natural log of the slope of the absorbance spectra between 275 and 295 nm by the natural log of the slope of the absorbance spectra between 350 and 400 nm. S R is a qualitative indicator of molecular weight of a sample, with larger values indicating lower molecular weight and smaller values having higher molecular weights (Helms et al., 2008). The specific UV absorbance at 254 nm (SUVA 254 ), which is the ratio of the decadic UV absorption at 254 nm divided by the DOC concentration, has been found to be useful in estimating DOM quality parameters such as: the aromaticity of DOM (Weishaar et al., 2003), the molecular weight (Chowdhury, 2013), and the source of DOM and type of environmental degradation (Hansen et al., 2016).

Dissolved Organic Carbon and Stable Carbon Isotope Analyses
For measurement of (DOC), samples were first sparged with ultrapure Argon gas for 20 min to remove dissolved inorganic carbon (DIC). Sparged samples were analyzed on an OI Analytical 1030D Aurora total organic carbon analyzer, using wet chemical oxidation, coupled to a Thermo Delta V Plus isotope ratio mass spectrometer (IRMS) to determine stable isotope values, expressed as the standard delta notation, δ 13 C-DOC, using Equation 2: where δ 13 C is the stable isotope ratio of carbon in parts per thousand, R sample is the ratio of 13 C to 12 C for the unknown sample, and R standard is the ratio of 13 C to 12 C for the known standard (Osburn and St-Jean, 2007).
(DOC) measurements were blank-corrected for ultra-pure Milli-Q water, then calculated using a linear regression curve of known caffeine standards with concentrations from 1 to 20 mg C L −1 . δ 13 C-DOC values were blank corrected and referenced to the Vienna Pee Dee Belemnite (VPDB) scale via a linear regression of six caffeine (IAEA-600, −27.7 ± 0.04 ) and two sucrose (IAEA-C6, −10.8 ± 0.03 ) International Atomic Energy Agency (IAEA) standards. Precision for (DOC) and δ 13 C-DOC values were ± 0.4 based on reproducibility and calibration to sucrose standard. Milli-Q blanks were run every 10 samples for quality control. DOC stable isotope values for different sources of natural OM from previous studies were tabulated for reference (Supplementary Table S3).

Geospatial Wetland Flood Model
The area and type of freshwater wetlands that were flooded in response to Hurricane Matthew were determined using Esri's ArcMap 10.5.1. Input variables included: (1) a high-resolution Light Detection and Ranging (LiDAR) digital elevation model (DEM) raster dataset with 6.3 cm vertical accuracy and 1 m horizontal accuracy downloaded from the NOAA Digital Coast data viewer (OCM Partners, 2017), (2) USGS gauge height and discharge rates, and (3) a coastal wetland GIS polygon vector dataset (Sutter, 1999). The relative elevation model (REM) tool within the Riparian Topography Toolbox (Wall, pers. comm.; Dilts et al., 2010) was used for removal of the slope trend of the floodplain from the DEM. Removal of slope trend of the floodplain from the REM was required because the slight slope of the ground to coastal water surface (1:1.8 m over the area of our survey) would cause the lower wetlands to become more flooded than the upper wetlands. Once removed, the water surface at base flow in the river and upper estuary was set to 0 m (REM) to flood an otherwise linear surface.
The study area ( Figure 1C) was extracted from the DEM using the Spatial Analyst toolbox (ESRI, 2011). The extracted wetland DEM was also processed by the REM toolbox to remove the slope trend from the surface water using kernel density with a search radius set to 3000 meters; this search radius was tested and found to be dependent on the size of the floodplain. The 3D Analyst toolbox (ESRI, 2011) was used to determine six surfaces of the flooded area and volume of the detrended freshwater riparian zone during Hurricane Matthew: relative 0 m, National Weather Service (NWS) flood stage, moderate flood stage, major flood stage, and maximum gauge height. A surface was also determined for mean gauge height over the period of October-December 2016, after the storm. Gauge height and discharge rates were obtained from the USGS gauging station on the Neuse River near Fort Barnwell, NC, United States (#02091814) (Supplementary Table S4). The value from the DEM at Ft. Barnwell (1.770 m) was removed from the gauge height values to account for removal of slope trend in the REM. In ArcMap, a constant planer raster layer was created (Spatial Analyst toolbox; ESRI, 2011) using the gauge height value, then the Cut Fill tool (3D Analyst toolbox; ESRI, 2011) was used to subtract the constant layer from the DEM. The Cut Fill raster layer was then extracted by attributes (Spatial Analyst toolbox; ESRI, 2011) where the volume was negative, to obtain the volume (V T ) and area (A T ) of the total flooded surface. An ArcMap model builder workflow plot indicates the processes that were used to determine the flood model (Supplementary Figure S1).
Next, North Carolina's Department of Coastal Management spatial wetland mapping layer (Sutter, 1999) was overlain on the flooded surface and extracted by mask (Spatial Analyst toolbox; ESRI, 2011) to determine the area (A W ) and volume (V W ) of wetlands present during various flood stages. The cell size of the DEM was 1.5 × 1.5 m. A W was calculated by multiplying the cell area (2.25 m 2 ) by the cell count of the total wetland flooded area from the attribute table in ArcMap. DOC mass loading (L) to the NRE from the Neuse River Watershed was calculated from the volumes determined in the ArcGIS flood model (Supplementary Figure S1) by first subtracting the V W from the V T (Eq. 3) to obtain the volume of the non-wetland area (V nW ), then the DOC mass loading from the wetlands (L DOCW ) was determined by multiplying the flooded volume of the wetlands by the DOC concentration from the wetlands (C W ) and dividing by 1000 to get units of kg C (Eq. 4). Next, riverine DOC mass loading (L DOCR ) was calculated by multiplying the DOC concentration from Fort Barnwell, NC, United States (riverine site; C FB ) by the non-wetland flooded volume and dividing by 1000 to value in units of kg C (Eq. 5). Finally, total DOC mass loading (L DOCT ) was determined by adding the DOC mass loads of the wetlands and the river and then converting to Gg C (eq. 6).

Bayesian Monte Carlo Mixing Model
A Bayesian Monte Carlo (BMC) three end-member mixing model was used to validate the geospatial model created during this study (Arendt et al., 2015). BMC has been shown to be useful in determining fractional contributions of unique endmembers to bulk samples in many different Earth surface systems, including glacial meltwaters (Bhatia et al., 2011;Arendt et al., 2015), soil nutrients (Chadwick et al., 1999), and stable and radiogenic isotopes in seawater (Pichler, 2005;Rickli et al., 2010). Within our estuarine and sound system, we identified the possible contributing end-member sources as wetland, river, and ocean water. We assumed that there were no other major end-members for this study. Therefore, the sum of the fractional contribution (f) of wetland (w), river (r), and ocean (o) for our estuary and sound samples is given in Eq. 7.
To determine the relative contribution of our unique end members (w, r, or o) to the bulk water (NRE or PS), we initially identified both SUVA 254 and δ 13 C-DOC as geochemical parameters with ranges specific to our end-members. However, due to linearity between SUVA 254 and δ 13 C-DOC values (Osburn et al., 2019a), a third component (S R ) was used to properly constrain the model. Thus, the combination of SUVA 254 , δ 13 C-DOC, and S R geochemical components allow us to define unique end-member geochemical compositions based on measurements from wetland (n = 22), river (n = 4), and ocean (n = 2; values from Atar, 2017) samples. End-member compositions are defined as mean SUVA 254 , δ 13 C-DOC, and S R values ± the standard deviation for each end member: where the wetland endmember composition is −29.41 ± 0.73 , 4.94 ± 0.50 L mg C −1 m −1 , and 0.83 ± 0.06; the river end-member composition is −25.77 ± 0.10 , 3.64 ± 0.61 L mg C −1 m −1 , and 0.88 ± 0.06; and the ocean end-member composition is −23. 01 ± 0.50 , 1.40 ± 0.50 L mg C −1 m −1 , and 1.82 ± 0.50, respectively. The incorporation of standard deviations within each end-member composition allows the model to account for natural variability and uncertainty associated with these end-member compositions.
Once the end members were geochemically defined, the BMC model was fit to our raw sample SUVA 254 , δ 13 C-DOC, and S R data from Neuse River Estuary (NRE; n = 70) and Pamlico Sound (PS; n = 52) samples were input to our BMC model to ascertain the f w , f r , and f o for each sample. Intuitively, the NRE and PS samples differ in SUVA 254 , δ 13 C-DOC, and S R compositions, with mean values ± standard deviations of −28.56 ± 1.13 , 4.12 ± 0.53 L mg C −1 m −1 , and 0.84 ± 0.06 for all NRE samples, and −27.46 ± 0.89 , 3.62 ± 0.32 L mg C −1 m −1 , and 0.99 ± 0.08 for all PS samples. The BMC model is able to provide further insight to these differences and calculate the most likely relative fractional contribution of each end member (f w , f r , f o ) to the NRE and PS samples based on the unique 3-component end-member geochemical compositions and associated uncertainties we defined. The BMC model uses prior and posterior probability density functions based on the known (end-member compositions and bulk NRE and PS compositions) and the unknown (relative contributions from each end member to the bulk samples) to determine the likelihood of any fractional contribution outcome. Each outcome is accompanied by a known uncertainty that places constraints on how likely each outcome is. The likelihood of a model prediction was calculated by Eq. 8.
Where subscripts p and o are the predicted and observed measurements of each variable, while σ is the measurement of uncertainty. One hundred million prior and posterior samples were tested for each bulk water sample and three standard deviations of the data for the bulk samples were used. The model was run separately for the NRE and the PS using the same defined end-member ranges. Models were run separately due to spatial and source differences between the NRE and PS, with NRE generally being more terrestrial sourced and PS more marine sourced. Average acceptance rates of the mixing model were determined via the ratio of prior samples to accepted posterior samples (see Eq. 8) and these rates were used to ascertain our confidence in the model fit.

S R
Mean daily S R values in the NRE increased from a low value of 0.77 ± 0.02 (17 October) to a maximum mean daily value of 0.90 ± 0.05 on 28 November, before decreasing slightly to 0.86 ± 0.09 on 13 December (Supplementary Figure S2C). NRE mean daily S R values were initially smaller than the mean river and wetland values (0.88 ± 0.06 and 0.83 ± 0.06, respectively), however, as the discharge began to decrease and wetland connectivity decreased the S R values became larger than the river and wetland values (Supplementary Figure S2C)

Flooding of Wetlands in the Lower Neuse River Watershed in Response to Hurricane Matthew
Hurricane Matthew's rainfall was localized well upstream of the NRE-PS, yet caused major flooding in the river's lower watershed above head-of-tide, its estuary, and the Pamlico Sound. Heavy precipitation from Hurricane Matthew in the Neuse River watershed occurred over a 2-day period from 08 October to 10 October (e.g., maximum of 419 mm at Kinston, NC, United States) (Figure 1A; Stewart, 2017). The precipitation associated with the storm caused extensive flooding throughout the watershed (Musser et al., 2017;Stewart, 2017). Flooding was most prevalent in the mid to lower Neuse River watershed, with flood recurrence intervals ranging from ca. 100 to 500 years (Musser et al., 2017). The USGS gauging station on the Neuse River at Fort Barnwell recorded gauge heights above the US National Weather Service Flooding of semi-disconnected to disconnected freshwater wetlands adjacent to the Neuse River's main channel between Fort Barnwell and ModMon sampling station NR0 at Streets Ferry Bridge was evident in satellite imagery comparing May and October 2016 (Figure 3). At low gauge height and low discharge, the wetlands were inundated with water; however, they were not connected via surface water to the Neuse River ( Figure 3A). Groundwater connectivity was assumed, but not known, and is not considered in this study. During the high gauge height and discharge associated with Matthew's rainfall, these wetlands became hydrologically connected to the Neuse River, which allowed for flushing of the water and OM present in these wetlands downstream into the NRE-PS (Figure 3B). This connectivity was apparent when examining results of the REM flood model (Figures 4, 5). During base flow and low gauge height conditions (0 m, REM) some wetlands were inundated (brown areas); however, the river was not connected via surface flow to these areas. The total area and volume that was flooded during base flow was 1.88 × 10 7 m 2 and 2.11 × 10 6 m 3 , respectively, while the flooded wetland area and volume were 1.55 × 10 7 m 2 and 1.76 × 10 6 m 3 . As the gauge height at Fort Barnwell, NC, United States increased, the areal and volumetric extent of wetlands that were inundated also increased. NWS flood stage (3.598 and 1.828 m; actual and relative gauge height, respectively) is the gauge height at which hydrologic connectivity of the wetlands to the Neuse River occurred (Supplementary Figure S3). A T and V T were 2.7 and 33.8 times larger during NWS flood stage than at base flow (5.03 × 10 7 m 2 and 7.12 × 10 7 m 3 , respectively). Similarly, increases in the A W and V W were present during NWS flood stage (2.6 and 35.3 times larger than base flow, respectively). A T , V T , A W , and V W all increased similarly in the moderate and major flood stage models as they did in the flood stage model ( Table 1). The V T and V W for the major flood stage were two orders of magnitude higher than at relative base flow ( Table 1). Moderate and major flood stage models also had more areas of inundated ground surface that was not classified as wetlands (Supplementary Figures S4, S5).
Following the passage of Hurricane Matthew over eastern NC, the peak of the flood pulse on 16 October 2016 caused the most flooding of the wetlands. The A T was 4.4 times greater than base flow, while the V T was two orders of magnitude higher than base flow (Table 1 and Figures 4, 5). Likewise, the A W (4.65 × 10 7 m 2 ) was ca. 3 times larger than base flow, and the V W (1.57 × 10 8 m 3 ) was two orders of magnitude greater than base flow (1.55 × 10 7 m 2 and 1.76 × 10 6 m 3 , respectively). Also, the area of flooded land that was not considered wetlands was greatest at this gauge height ( Figure 5). Comparatively, the A T , V T , A w , and V W at the mean gauge height from 01 October to 31 December (1.930 m, actual; 0.160 m, REM; Supplementary Figure S6) was more similar to base flow than during the flood pulse ( Table 1).

Model Estimates of DOC Fluxes Due to Hurricane Matthew
Flooding following the passage of Hurricane Matthew caused hydrologic connectivity of the riparian wetlands in the Neuse River's floodplain with the main channel, and substantial export of DOC from the riparian wetlands into the lower Neuse River downstream and into the NRE and PS.
Once volumes were determined for the A T and the A W , the estimated DOC fluxes from the riparian wetlands and from the river into the NRE and PS were calculated using Eq. 3-6 ( Table 2). Carbon mass loading from the wetlands was assumed to occur from the start of wetland connectivity until the wetlands became disconnected after the flood pulse had passed. This period was ca. 31 days (08 October through 08 November 2016). Connectivity was determined to occur at a gauge height above 0.7 m on the REM (actual gauge height 2.47 m). The gauge height at each day was matched with the model from different flood stages (2-days, maximum flood stage; 5-days, major flood stage; 3days, moderate flood stage; 3-days, flood stage; 18-days, 0.7 m REM). These values were summed to determine the total DOC loading to the NRE.
The L DOCR was determined to be 3.76 Gg C, while the L DOCW was 12.65 Gg C, thus the total L DOCT was 16.41 Gg C ( Table 3). Osburn et al. (2019a) determined the DOC stock change from the storm in the estuary using only concentrations and volume estimates was 4.74 Gg C. Therefore, the modeled L DOCT was ca. 12 Gg greater than the previously calculated DOC stock change. However, this discrepancy in DOC mass loading can be explained by the calculations for the DOC stock change from Osburn et al. (2019a) only covering 1 week of the sample period, while the L DOCT computed here using the flood model was done over the entire connectivity period (31-days), once that duration was known. Multiplying the DOC stock change in the NRE by 4 gives a value comparable to our modeled number (18.96 vs. 16.41 Gg C, respectively). For context with annual loads in this system, the hydrologic load estimator model, LOADEST (Runkel et al., 2004) was used to determine the estimated annual DOC load from the Neuse River during 2016 (51.32 Gg C) by using discharge data from Fort Barnwell combined with DOC concentrations throughout the year. Combined, L DOCR and L DOCW were a large proportion (ca. 32%) of the annual estimated load (7.3 and 24.6%, respectively; Table 3). It was clear from these model calculations that wetlands supplied the bulk of the DOC exported to the NRE-PS as a result of Hurricane Matthew.

Bayesian Monte Carlo Mixing Model of DOM Sources
The NRE was assessed both spatially and temporally using the BMC mixing model. The acceptance rates for our three endmember BMC mixing models ranged from 0.15 to 0.48 for the NRE and 0.31 to 0.49 for the PS. The average acceptance rates of each model were 0.30 and 0.41 (NRE and PS, respectively), which indicated that the PS model performed better than the NRE model ( Table 4). Rosenthal (2011) suggested that BMC mixing models would provide accurate results if the average acceptance rates were between 0.1 and 0.6, while Lunn et al. (2009) indicated that average acceptance rates should be between 0.2 and 0.4, with 0.4 being optimal.
Similar to results from our geospatial model, the BMC geochemical mixing model showed that the wetlands were the largest contributor of DOC to the NRE and PS following Hurricane Matthew. Spatially across the entire NRE, the mixing model indicated the f w of DOC was the dominant contributor at ca. 40-90%, followed by the f r (10-50%), with the ocean end-member generally contributing the least ( Figure 6A). Temporally, throughout the period of wetland hydrologic connectivity (08 October through 08 November), the source of DOC to the NRE was wetland dominated (Figures 2, 6). Once the river returned to baseflow in late November, there was a transition to a more mixed water source with the river source dominating the DOC (Figures 2, 6).
Examining the source fraction contributions to the estuary spatially showed how the wetland source decreased downstream (Figure 6). The estuary sites were spatially divided into upper, middle and lower sections. The upper section included sites NR0-NR50, while the middle section included sites NR60-120 just before the bend in the estuary (Figure 1B). The lower section included the remaining sites downstream of the bend in the estuary into southern Pamlico Sound (NR140-NR180). The f w was predominant in each section of the estuary with few exceptions. In the upper section, there were two sample dates where the river and wetland contribution indicated mixing (28 November and 13 December), after hydrologic connectivity occurred ( Figure 6B). For all spatial sections of the estuary, the f r and f w overlapped on 28 November during the lowest flow of the sampling period (Figures 2, 6B-D).
The fractional contribution of wetlands, river, and ocean for PS sites were only examined temporally, due to the spatial proximity of the PS sampling sites (Figure 1B). Within the PS, the f w ranged from ca. 35-65%, which was lower than fractions in the estuary, due to continued dilution of this source from the upstream wetlands. During the initial sample day (27 October) ca. 3 weeks after the storm, the f r was the largest contributor to the PS (Figure 7). However, once the initial flood pulse reached the PS (02 November) the wetland became the larger contributor. On 09 November mixing of these three sources was apparent due to overlap between the f r , f w , and f o . The remaining sampling period indicated that wetland DOM was dominant in the PS (Figure 7).
Because the BMC mixing model accounts for natural variance in end-member compositions, measurement uncertainty, and factors in the variance and covariance of geochemical indicators, it produces uncertainties associated with each most likely fraction contribution of each end member. The uncertainties associated with the estimated fractional contributions are represented as error bars in Figures 6, 7 and decrease with the degree to which the end-member geochemical compositions are distinguishable from one another. For NRE and PS samples with overlapping end-member fraction contribution estimations, Figures 6, 7 still show the most-likely fraction contribution but overlapping error bars indicate there is a much less-likely possibility that the fractions exist within the overlapping space. However, the observed fraction contribution trends for our NRE and PS samples and findings from our geospatial model strongly indicate the BMC most likely contributions are likely valid.

DISCUSSION
The main hypothesis of this study was that flooding of freshwater riparian wetlands surrounding the Neuse River, would become a major component of DOM in the NRE-PS due to hydrological connectivity. This hypothesis was confirmed in two ways. First, by expanding on an existing flood model to quantify the areal  extent and volume of flooded wetlands via their hydrological connectivity to the main stem of the river, we determined the volume of water (and then using concentrations, DOC) exported from the wetlands to the NRE-PS (Table 4 and Figures 4, 5,  Supplementary Figures S3-S6). Second, this study combined stable C isotope, SUVA 254 , and S R values in a Bayesian Monte  Annual load of DOC was estimated using Load Estimator (LOADEST) software (Runkel et al., 2004). Percent of annual load was calculated by dividing the estimated annual load by the calculated load. DOC mass flux from the terrestrial environment was determined by dividing the DOC load by the number of days the wetlands were hydrologically connected. n.a.-not applicable. Carlo mixing model to determine estimates of the wetlands' contribution to DOM in both the NRE and PS (Figures 6, 7). This study revealed the substantial contribution of riparian wetlands to the DOC flux from a coastal watershed into receiving waters following a major hurricane. Within the context of the Pulse-Shunt concept, we discuss how our new combination of geospatial and geochemical modeling provides solid estimates on not only the magnitude of DOC fluxes from extreme events but also insight into where the DOC originates. Combined, these modeling approaches offer a key means to track source and fate of DOC delivered to coastal waters after an extreme event. Implications of the exported DOC in coastal waters are discussed.

Pulse: Wetlands as a Primary Source of DOC to Coastal Waters Following an Extreme Weather Event
Applying the PSC to our results, the change in DOM quality in the estuary and sound was attributable to the pulse of wetland DOM caused by their flooding and hydrological-connectivity to the Neuse River. Using the flood model, it was determined that initial connectivity of the wetlands to the main channel of the Neuse River occurs between a gauge height of 2.274-2.774 m (0.5-1.0 m, REM), though the exact level is variable between the upper and lower wetlands. For context, mean annual gauge height at Fort Barnwell has ranged from 1.32 to 1.98 m over the past 17 years, while mean annual discharge at Fort Barnwell has ranged from 50 m 3 s −1 to 181 m 3 s −1 over the past 21 years (U.S. Geological Survey [USGS], 2018). Therefore, the increase in downstream export of wetland OM could occur prior to NWS flood stage at Fort Barnwell (3.598 m). This metric is important to establish for coastal rivers with respect to wetland pulses in terrestrial DOC export. Spencer et al. (2013) showed very strong correlation (R 2 = 0.81) between export of CDOM as a marker for terrestrial DOC, and the percent of wetland in a watershed. Hanley et al. (2013) found high positive correlation (R 2 = 0.91) of watershed %wetland cover with SUVA 254 values in coastal rivers. The watershed% wetland in the Neuse River watershed is 13% (Osburn et al., 2016), which is comparable to the wetlandrich eastern North American rivers (e.g., the Edisto and the Androscoggin) which export substantial amounts of CDOM to the ocean (Spencer et al., 2013).
Hence our results strongly suggest the influence of wetlands in the DOM characteristics of the Neuse River, its estuary, and Pamlico Sound. While the range of our riverine δ 13 C-DOC values (−25.9 ± 0.6 ) were similar to results reported from other rivers in the literature including the Tar (NC) and Santa Clara (CA) Rivers (Supplementary Table S3), the relatively depleted estuarine δ 13 C-DOC values and high SUVA 254 values were consistent with observations of the riparian wetlands. This shift to depleted isotope values is consistent with flooding: in the lower Mississippi River, where for example, δ 13 C-DOC values decreased from −25.5 to −27 as discharge increased (Cai et al., 2015).
Finally, comparing the flood and geochemical modeling approaches revealed similar estimates of wetland contributions to DOC in the estuary and sound, therefore providing independent measures of the pulse. The flood model estimated that 77% of the DOC load from Matthew was due to wetlands ( Table 3). This model was purely physical in the sense that DOC concentration and water volumes were the two parameters used to compute relative fractions of wetland and river in the DOC load. The BMC mixing model estimated that 70 ± 20% of NRE DOM and 52 ± 20% of Pamlico Sound DOM exhibited a wetland character, considering only mixing of river, wetland, and oceanic sources in both systems ( Table 4). Segmentation of the NRE into upper, middle, and lower regions further attested to the strong influence exerted on DOM properties by wetlands (Figures 6B-D). Segmentation also provided a view of the time-lag of the wetland influence on DOM in the estuary. For example, wetland contribution in the middle estuary peaked at ca. 85% on 08 November, falling to roughly 50% 2 weeks later ( Figure 6C). In the Pamlico Sound the% wetland contribution was consistent at ca. 60% from 02 November through 14 December ( Figure 6D). The segmentation and time-series results indicated that the wetland characteristics remained a prominent feature of the NRE and PS for nearly 3 months following Matthew (Osburn et al., 2019a).

Shunt: Rapid Delivery of a Majority of DOC Coastal Waters in Response to an EWE
Inland precipitation related to Matthew resulted in a "shunting" of the wetland "pulse" of OM into the NRE-PS which dramatically increased its stock of DOC. Use of long-term monitoring data sets such as ModMon in the NRE-PS allowed for comparisons to be made between extreme events and long-term averages in these systems. We combined analyses from this study with long-term trends in DOC to determine what effect Hurricane Matthew's floodwaters had the shunt of DOM to the NRE-PS. Osburn et al. (2019a) indicated there was ca. 65% increase in DOC stock in the NRE-PS in response to Hurricane Matthew. When comparing our results from the geospatial and geochemical models, we observed that both produced similar fractional estimates for the river (23%) and wetland endmembers in the NRE (70 and 77% for the geochemical and flood models, respectively). The similarity in these estimates is encouraging and demonstrates the need to incorporate lateral flooding of wetlands in coastal watersheds to provide best estimates of C export from land to coastal waters.
About 30-70% of annual DOC export from watersheds occurs following storms (Raymond and Saiers, 2010). Hurricane Matthew's floodwaters caused a total terrestrial DOC load of 16.41 Gg C between 08 October and 08 November, during the period that the riparian wetlands were hydrologically connected to the Neuse River. This episodic load of DOC accounted for ca. 32% of the estimated annual load of DOC for the NRE (Table 4) which is similar to what other studies have shown in larger systems. For example, flooding from storm events in the Lower Mississippi River Basin in 2011 caused the DOC export from the Atchafalaya River to coastal waters to increase from ca. 5 to ca. 22 Gg C d −1 . This increase accounted for ca. 50% of the annual carbon load . Our results indicate that due to the hydrologic connectivity of the wetlands upstream of the estuary, around 70-77% of the 32% annual load increase in the NRE can be attributed to the wetlands. Similarly, Majidzadeh et al. (2017) indicated that in the Yadkin-Pee Dee Hurricane Matthew exported a large amount of wetland DOC.
Our flood model estimates of DOC flux were conducted using a different approach than prior work, which strongly suggests a convergence on accurately estimating the DOC load caused by storms. For example, Paerl et al. (2018) Hirsch et al., 2010;Hirsch and De Cicco, 2015) model. When compared to other storm events in the NRE the present DOC mass loading results were higher than estimated DOC loads from Hurricanes Fran (1996) and Irene (2011) although much lower than estimates from the combination of sequential Hurricanes Dennis, Floyd, and Irene in 1999 (7.4-8.4, 14 (TOC), and ca. 60 Gg C, respectively; Brown et al., 2014;Crosswell et al., 2014;Paerl et al., 1998;Paerl et al., 2018).
The Pulse-Shunt concept proposed by Raymond et al. (2016), was based on headwater fluvial systems but describes the dynamic we have modeled for coastal rivers and estuaries. First, during high discharge, from a low-frequency, large storm event, regional flooding inundates wetlands, with the hydrologic-connectivity of wetlands to the main channel of the river facilitating a flushing of DOM into the river and downstream to receiving waters. Second, as for lower order streams and rivers, estuaries could be turned into passive pipes. Evidence for the NRE existing in a passive pipe mode comes from characteristic where up to ca. 40% of the river's annual load of terrestrial DOM was shunted downstream into Pamlico Sound. Moreover, this last point is important; this study demonstrated the lingering effect of these episodic and infrequent storms and suggest that coastal ecosystems such as the NRE attenuate the floodwater pulse over long periods of time.

Process: Significance of Extreme Weather Events on Coastal Carbon Cycles
Our findings demonstrated that EWE such as Hurricane Matthew, the second 500-year storm to occur in NC within the past 20 years , transfer substantial amounts of terrestrial C to coastal waters in just a few weeks. Osburn et al. (2019a) provided evidence suggesting that much of this material can be processed and sustain Pamlico Sound as a weak CO 2 source to the atmosphere. Increases in stock of relatively "fresh", sunlight-absorbing DOC into coastal waters increases the potential for it to be converted to atmospheric CO 2 via microbial respiration and photooxidation during a long residence time (Spencer et al., 2009;Bauer et al., 2013;Ward et al., 2017). Conversion of allochthonous DOM that would normally be stored in the terrestrial environment in soils or wetland sediments to CO 2 in coastal waters creates a feedback loop to the atmosphere, which potentially could influence coastal C cycling (Bauer and Bianchi, 2011).
The geomorphological properties at the land-ocean continuum will modulate the transport and fate of DOM in coastal waters (Bianchi et al., 2007). For example, the NRE-PS is an embayed estuarine system which allows for some export of OM to the coastal Atlantic Ocean through three small inlets. However, the Mississippi and Amazon rivers are open deltas that permit OM to be transferred out into the open ocean (Hedges et al., 1997;Bianchi et al., 2013). Residence time in large river plumes may approximate the long residence time (ca. 1 year) within Pamlico Sound (Supplementary Table S1). These differences suggest that in similar embayed lagoons and semi-enclosed estuaries the fate of a large portion of the OM that would normally be stored in the estuaries or exported to the coastal ocean (e.g., DOM in large river plumes) due to flooding caused by extreme precipitation events, is alternatively processed in situ. Long residence time lagoonal estuaries constitute about 25% of estuaries globally and are a predominant feature of the SE Atlantic and Gulf coasts of the United States (Dürr et al., 2011).
Extreme weather events are increasing in intensity due to a positive feedback loop created by global warming that will undoubtedly cause a redistribution of water across the planet (Trenberth, 2011). Increased heating leads to longer periods of drought (Trenberth, 2011); extensive dry periods also allows terrestrial OM to accumulate in hydrologically isolated wetlands. However, warming also increases atmospheric water vapor, which leads to a larger moisture supply to fuel more intense precipitation events (Trenberth, 2011) and thus more flooding events . Large pulses of reactive terrestrial OM to coastal waters have the potential to contribute to the increase in atmospheric CO 2 through photooxidation and microbial respiration of DOC in estuaries (Bauer et al., 2013). Saturated soils from heavy rainfall in the Neuse River Watershed prior to Hurricanes Fran and Matthew caused similar exports of total organic carbon to coastal waters (ca. 14 and 21 Gg C, respectively), indicating the importance of extreme precipitation events and flooding on the transference of terrestrial OM to the coastal ocean. We suggest that these weather extremes (dry periods punctuated by intense precipitation events) will amplify pulses of C stored in wetlands which then are shunted into receiving coastal waters where this organic C is returned to the atmosphere as CO 2 . A predicted rise in the frequency of EWE, including major hurricanes in the Atlantic, where many lagoonal estuaries exist, will continue to exacerbate these conditions (Bender et al., 2010;Trenberth, 2011;Lehmann et al., 2015;Walsh et al., 2016).

CONCLUSION
A Bayesian Monte Carlo (BMC) mixing model combining stable carbon isotope analysis, SUVA 254 values and S R values, revealed new insights into the effects of extreme events on this system, improving upon prior work in which OM sources from events such as hurricanes had not been clearly identified (Peierls et al., 2003;Osburn et al., 2012;Brown et al., 2014;Paerl et al., 2018). The BMC method provided estimates of the fractional source contributions from wetlands to the NRE following Hurricane Matthew in 2016, which we then compared to estimates made with a hydrologicalconnectivity flood model. Through this comparative modeling approach, we determined that the riparian wetlands exported ca. 50-70% of the DOM to the NRE-PS during Hurricane Matthew's floodwaters, confirming predictions from earlier work (Osburn et al., 2019a).
Our modeling approach further provided insight into estuarine responses to low-frequency large storm events in a manner similar to rivers under the PSC (Raymond et al., 2016). For the NRE-PS, this coastal system responds similarly: the estuary essentially became a passive pipe, resulting in ∼40% of the annual load of DOM to be shunted through the estuary into Pamlico Sound. An important insight we gained from this approach was that the pulse of DOM from the hydrologically connected wetlands following Hurricane Matthew accounted for ca. 25% of the estimated annual load to the Neuse River Estuary. PSC could be amended for coastal environments such as Pamlico Sound to include processing (P) of DOM into CO 2 by photooxidation or microbial respiration. Our results suggest that EWE, such as tropical cyclones, have a significant effect on coastal carbon cycles, by inducing large pulses of terrestrial DOM from floodwaters associated with them and translocating OM accumulating over time in wetlands into biogeochemically dynamic ecosystems where substantial processing can occur (Crosswell et al., 2014;Osburn et al., 2019a).
Estuarine and coastal waters are home to more than half of the world's population and are important economic, recreational, and environmental ecosystems . However, a lack of high-resolution observational data hinders our comprehension of the effects of EWE on coastal biogeochemistry (Brown et al., 2014). Therefore, the development of new interdisciplinary methods such as what we have done for this study is crucial to understanding the coinciding and subsequent processes occurring in coastal waters in response to EWE.

DATA AVAILABILITY STATEMENT
The data used in this project can be found at the NSF Biological and Chemical Oceanographic Data Management Office (BCO-DMO; https://www.bco-dmo.org/project/734599).

AUTHOR CONTRIBUTIONS
JR participated in the experiment design, performed the field samplings of the river and wetlands, conducted laboratory analyses, performed the data analyses, and was responsible for the manuscript preparation. CA contributed to the analysis of the geochemical mixing model and participated in the manuscript preparation. AH participated in field samplings in the estuary and sound, conducted laboratory analyses, and participated in the manuscript preparation. HP provided the estuary and sound samples and participated in the manuscript preparation. CO participated in the experiment design, data analysis, and manuscript preparation.