Five Years Measuring the Muck: Evaluating Interannual Variability of Nutrient Loads From Tidal Flooding

Due to sea level rise, tidal flooding is now common in low-lying coastal systems around the world. Yet, the contribution of tidal flooding to non-point source nutrient loads and their impact on the quality of adjacent waters remains poorly constrained. Here, we quantified dissolved nutrient loading and Enterococcus abundance during annual autumnal king tides (i.e., perigean spring tides), between 2017 and 2021, in a sub-watershed of the lower Chesapeake Bay. To calculate nutrient loading from tidal flooding, we used geospatial inundation depths from a street-level hydrodynamic model to estimate floodwater volumes during each of the five sampling events and the difference between nutrient concentrations in floodwater and pre-flood measurements. Results showed that dissolved nutrient concentrations were higher in floodwaters than in estuarine waters and resulted in dissolved nitrogen and phosphorus loads that reached 58.4 × 103 kg and 14.4 × 103 kg, respectively. We compared our load estimates to the tributary-specific total and land-based federal allocations (i.e., total maximum daily loads (TMDL)) for total nitrogen (TN) and total phosphorus (TP). Even the more conservative calculations indicate that inputs of dissolved nutrients during a single tidal flooding event can exceed 100% of the annual load allocation. Additionally, more than 80% of the floodwater samples collected each year showed Enterococcus abundance that exceeded the threshold for recreational water use in Virginia (104 MPN 100 ml−1). Failing to account for non-point source loading of nutrients and contaminants from tidal flooding as sea level rises could result in worsening eutrophication and deterioration of coastal economies and the health of coastal communities around the world.


Introduction
Coastal ecosystems are under pressure from multiple humaninduced stressors including hydrologic and land use changes, eutrophication, and warming water and air temperatures (Day Jr et al. 2012;Lotze et al. 2006;Paerl et al. 2006a). In the USA and around the globe, one of the most significant problems in coastal waterways is the continued input of anthropogenic nutrients, including dissolved inorganic nitrogen and phosphorus from food production and processing (Pinckney et al. 2001;Smith 2003;Seitzinger et al. 2002;Scavia and Bricker 2006). These nutrients can be introduced via point source discharges and non-point sources including forested, agricultural, commercial, residential, and urban runoff (Sabo et al. 2022). Despite years of targeted nutrient reduction programs, primarily aimed at point source discharges, non-point source nutrient inputs remain difficult to control. These are largely due to increases in human populations in coastal areas and consequent increases in fertilizer use to support agriculture (Swaney et al. 2012) to feed those growing populations (Seitzinger et al. 2002). Excess nutrient inputs contribute to eutrophication, the increase in intensity and frequency of harmful algal blooms, and the expansion of "dead zones" where oxygen is depleted (Kemp et al. 2005;Rabalais et al. 2009).
The capacity for storms and storm-related flooding and runoff to transport pollutants to coastal estuaries and waterways has been recognized for decades (Magenis 1988;Hale et al. 2015;Wu et al. 2016;Paerl et al. 2006b), and the relationship between land use and non-point source nutrient loading has been examined for agricultural, urban, forested, and wetland landscapes (Basnyat et al. 1999;Paerl et al. 2006a, b;Swaney et al. 2012;Selbig 2016;Wu et al. 2016). However, in coastal areas, tidal flooding (even when the skies are blue) presents a potentially different mechanism for delivering nutrients to adjacent estuaries (Macías-Tapia et al. 2021;Weissman and Tully 2020). Unlike stormwater runoff, the flow is not unidirectional; rather, coastal waters advance and inundate the landscape during high tide and then retreat as the tide ebbs. The residence time of water on the landscape is longer, during which time floodwater can both mobilize and react with material accumulated on the landscape before returning to the estuary as the tide recedes. Additionally, previous studies have shown how interactions between floodwater and groundwater can modify the input of dissolved organic matter and inorganic nitrogen and phosphorus to the estuary (Wilson and Morris 2012;Santos et al. 2021). Given increases in the intensity, frequency, and magnitude of tidal flooding events over the last few decades due to sea level rise (Nicholls and Cazenave 2010;Spanger-Siegfried et al. 2014;Neumann et al. 2015;Ezer 2018), these inputs are likely to increase in the future as sea level continues to rise.
The impacts of coastal flooding on human health in urban areas have been well documented in many parts of the world (Hajat et al. 2005;Fernandez et al. 2015;Rui et al. 2018). However, as with nutrient loading, most studies are focused on storms and storm-induced flooding. Macías-Tapia et al. (2021) found that during a "blue sky" tidal flooding event in the Lafayette River watershed, a lower Chesapeake Bay sub-tributary, more than 90% of floodwater samples collected exceeded Virginia's recreational standard for Enterococcus abundance (104 MPN/100 ml; MPN indicates the most probable number), a fecal matter indicator. An increasing proportion of the population is interacting directly with flood waters as coastal flooding becomes more common, exacerbating human health concerns associated with floodwater contact (Neumann et al. 2015;Tellman et al. 2021). Furthermore, multiple studies have shown that flooding disproportionately affects low-income communities that are unable to relocate or pay for measures to reduce their flood risk (Maantay and Maroko 2009;Collins et al. 2018;Chakraborty et al. 2020).
Quantifying nutrient loading from tidal flooding is complicated. It requires knowing the inundation volume, considering the local hydrology, and knowing the difference between the nutrient concentrations in the estuary prior to flooding the landscape and collected floodwater samples. To develop models that can be applied at a watershed level, samples need to be representative of different land uses and flooding conditions. This requires multiple sampling events over many years. Here, we compare nutrient concentrations from floodwater samples collected during a perigean spring tide flooding event (a.k.a., king tides) held each autumn between 2017 and 2021. The goal of this study was to quantify the variability in dissolved nutrients and Enterococcus loading due to tidal flooding in the Lafayette River, a residential subwatershed of the lower Chesapeake Bay. If the magnitude of these loads is higher than recommended limits, they can exert significant threats to the health of local waterways and the health of those residing in impacted areas.

Methods
To quantify nutrient loads resulting from tidal flooding, it is necessary to quantify three parameters: nutrient concentrations in floodwaters, nutrient concentrations in estuarine waters prior to their encroachment onto the adjacent land, and floodwater volume. Because of the high variability of nutrient concentrations in floodwaters, we required a large number of samples to provide statistically robust estimates of these nutrient loads. So, to measure nutrient and bacterial concentrations in floodwaters, we coordinated annual citizen science sampling campaigns (named "Measure the Muck") during perigean spring tides (king tides) in each year from 2017 to 2021. These campaigns were enabled by a larger citizen science event (Catch the King Tide) organized by local media outlets in the Hampton Roads region of southeastern Virginia aimed at building public awareness about coastal flooding due to sea level rise. This event was made possible by the development of a phone application (Sea Level Rise) by a local non-profit group (Wetlands Watch) that provided a platform for citizens to collect data on water levels during tidal flooding events (Loftis et al. 2019).

Study Site
The Lafayette River is a sub-tributary of the lower Chesapeake Bay, and its watershed is largely residential and limited to the city of Norfolk, VA (Fig. 1). This city has one of the highest rates of relative sea level rise in the USA (Spanger-Siegfried et al. 2014), due to multiple factors including eustatic sea level rise, land subsidence, and changes in the strength of the Gulf Stream (Ezer and Corlett 2012;Ezer et al. 2013). This has resulted in a dramatic increase in the number of hours in which tidal flooding inundates the city each year over the last few decades (Ezer 2018).

Floodwater Sampling
The Measure the Muck (MTM) campaigns leveraged a citizen science program (Catch the King Tide) aimed at determining the extent of floodwaters during perigean spring tides annually from 2017 to 2021. Volunteers for the MTM sampling were recruited primarily from the Norfolk Public Schools and Old Dominion University (ODU). Prior to the Catch the King Tide event, volunteers were trained on how to use the Sea Level Rise phone application to record data and how to collect floodwater samples. On the day of the King Tide/MTM event, volunteers were armed with kits containing sampling supplies and sent to areas known to flood as described in Macias-Tapia et al. (2021). Once in a flooded area, volunteers mapped the extent of floodwaters using the phone application and collected metadata and samples near or just after high tide in floodwaters < 0.2 m deep. Detailed instructions were given during training sessions, and new or younger volunteers were accompanied by more experienced participants (e.g., staff or graduate students from ODU). While we had over 50 volunteers during the first three events, the number of volunteers dropped by more than half in 2020 and 2021 due to the COVID-19 pandemic. Despite differences in the number of volunteers, the number of samples collected was higher in the pandemic years (Table 1). This might had been due to the increasing experience on the organizers and volunteers. The fewest samples were collected in 2019 (n = 146), when flooding was limited in extent, and the greatest number of samples (n = 196) was collected in 2021 when the extent of flooding was greatest. During 2017During , 2018During , and 2020, and 180 samples were collected, respectively. Due to accessibility, most of the samples were collected along the southern shore of the estuary where there was more flooding on public property ( Fig. 1). For all the years, samples were collected from the mouth to the upper part of the Lafayette River.
To collect samples, the sampler donned reusable Playtex cleaning gloves and submerged the sample bottle (acid cleaned Nalgene™ polycarbonate bottles) into the floodwater to rinse and then fill. After collection, each sample was placed inside the provided coolers with ice packs. The sampling team recorded the sample number (each sample bottle issued was pre-labeled with a unique sample identifier), the GPS coordinates of the location the sample was collected, and a description of the sample site in both the Sea Level Rise application and on log sheets provided to the sampling team. After collection, samples were transported to a laboratory at ODU within 1 h of sample collection. Once in the laboratory, samples were processed for water quality analysis. Details on processing and analyses are provided below. Sample processing was supervised and accomplished by trained graduate students and university staff. However, as part of the public engagement component of the project, a subset of volunteers was trained to help with sample filtration, organizing samples, measuring salinity with a refractometer, and labeling sample tubes. 1 3 During each MTM campaign (i.e., 2017 to 2021), 40 samples were collected to determine Enterococcus abundance in floodwater. Samples were collected into sterile glass bottles, placed in a cooler with ice, and returned to the laboratory. Sample analyses were initiated within 24 h of their collection.

Estuarine Sampling
Prior to flooding events, samples were collected from the Lafayette River so that floodwater nutrient concentrations could be corrected for ambient nutrient concentrations. For the rest of the document, the term "baseline" is used to refer to conditions in the estuary prior to tidal flooding. In order to capture the spatial and temporal variability of the biochemical parameters previous to each flooding event, three types of baseline data were used in this study ( Table 2). The temporal and spatial differences among baseline datasets are expected to create different responses when compared with values from floodwater samples. The four baselines were (1) data from a time series site in the Lafayette River (time series), (2) long-term data from a Chesapeake Bay Program monthly water quality monitoring station in the Lafayette River (CBP), (3) CBP data collected just prior to each flooding event (CBP-year), and (4) data collected during cruise transects along the Lafayette River mainstem made in the days prior to the onset of tidal flooding (transect).
For the time series baseline, between 2017 and 2021, nutrient concentrations at the Lafayette River time series station were measured almost daily throughout the summer and early fall (Table 2). Samples were collected using a custom-made pump profiling system comprised of a Masterflex L/S 25 peristaltic pump and platinum-cured silicone tubing marked at 25-cm intervals with electrical tape. The submerged end of the tubing had a weight attached to keep it vertical while sampling. At the sampling end of the tubing, a 0.2-µm Whatman polycap filter was inserted into the sample stream to collect filtrate to measure dissolved constituents. Samples were filtered directly into sterile Falcon® tubes, placed into a cooler with ice packs for transport to the laboratory (< 1 h), and stored frozen (− 4 °C) until analysis (see below). Near-surface samples (< 1 m depth) were collected in this way at two stations in the Lafayette River ( Fig. 1), one near the headwaters and the other near the mouth. While samples were collected throughout the summer at these stations, only data from late summer and fall were used to provide in-estuary nutrient estimates for this study. The number of samples available per year varied with 117, 51, 102, 27, and 39, from 2017 to 2021, respectively (Table 2).
Long-term data for nutrient concentrations in the Lafayette River are available through the Chesapeake Bay Program. These data are publicly available (Chesapeake Bay Program Water Quality Data Hub 2021). Nutrient samples have been collected at stations near the mouth and head of the estuary since 1988 (Fig. 1). Sampling frequency has varied over the years at these sites. In this study, we used data from 1998 to 2018 when data collection was consistent at both stations within the Lafayette River (Table 2), and we used only concentrations measured in surface water samples collected above the pycnocline. From this dataset, we produced climatologies for DIN, DON, and PO 4 3− for the months of October and November, the months in which the king tide sampling events were conducted between 2017 and 2021. This data is hereafter referred to as Chesapeake Bay Program (CBP) data. Because nutrient concentrations can vary over short timescales in this system (e.g., Morse et al. 2014), in addition to comparing floodwater and average baseline concentrations, we made comparisons between floodwater concentrations and the most closely paired in-estuary nutrient concentration measurement made as part of the Chesapeake Bay Program monitoring for each year of Measure the Muck (CBP-year). The number of values available from CBP and CBP-year was 87 and 1 for all years, respectively (Table 2).

3
Samples for the boat transect baseline were collected during 2019, 2020, and 2021, from cruise transects from the mouth to the head of the Lafayette River (Fig. 1). The number of samples collected varied each year from 27 in 2021 to 84 in 2019 and 2020 with n = 39 (Table 2). Samples were collected several days before the peak of the perigean spring tide prior to the onset of flooding (Fig. 1). An identical pump profiling system as that used at the time series site was used to collect near-surface water samples while on the boat. As for the time series site sampling, a 0.2-µm Whatman polycap filter was affixed to the sampling end of the tubing, and nutrient samples were filtered directly into sterile Falcon® tubes, stored in a cooler with ice, and transported to the laboratory where they were frozen until analysis.

Sample Processing and Analysis
Once floodwater samples were returned to the laboratory, they were filtered under low vacuum pressure (< 5 in Hg) through Advantec GF75 filters (pore size of ~ 0.3 µm). The filtrate was collected into sterile Falcon® tubes and stored in a − 4 °C freezer until analysis. Dissolved ammonium (NH 4 + ), nitrite (NO 2 − ), nitrite + nitrate (NO x ), total dissolved nitrogen (TDN), and phosphate (PO 4 3− ) concentrations were measured as described in Macias-Tapia et al. (2021). Briefly, NO 2 − , NO x , and PO 4 3− were measured using colorimetric methods on an Astoria Pacific nutrient autoanalyzer according to the manufacturer's specifications (Hansen and Koroleff 1999). Ammonium concentrations were measured using the phenol hypochlorite method (Solorzano 1969). To measure TDN concentrations, samples were pre-treated with potassium persulfate to oxidize all of the dissolved nitrogen to nitrate and then measured as described above. The oxidation was done during a 30-min liquid cycle in an autoclave at 120 °C. Dissolved inorganic N (DIN) was calculated as the sum of NH 4 + and NO x concentrations. Dissolved organic N (DON) concentrations were calculated as the difference between TDN and DIN concentrations. For each analyte, the detection limit (DL) was calculated as three times the standard deviation of the lowest point of the standard curve (i.e., 3*σ, n > 3). We assigned the calculated DL to samples that showed concentrations lower than that. The same chemical and numerical analyses were used to measure dissolved nutrient concentrations in floodwater and estuarine baseline samples.
Enterococcus abundance was measured using Enterol-ert® kits (IDEXX Laboratories, Inc.) following the manufacturer's specifications. Briefly, this method calculates the most probable number (MPN) of Enterococci bacteria based on a statistical probability analysis of the number of blocks that activate the fluorescent indicator (Ramoutar 2020). Blanks of autoclaved pure water were used to account for contamination during analysis. Results for Enterococcus abundance in floodwater were compared with standards established by the US Environmental Protection Agency (EPA) and the Virginia Department of Health (VDH) for safe recreational use. Virginia adopted the federal standard of 104 MPN ml −1 as the limit for Enterococcus abundance in natural waters deemed safe for recreational use (EPA 2010).

Statistical Analysis
For both floodwater and baseline datasets, we calculated the arithmetic means and standard deviations (SD) and the medians for each parameter. Before performing any statistical comparisons between floodwater and baseline concentrations of analytes, we used the Kolmogorov-Smirnov test to check for normality in the different datasets. For most of the datasets (excluding 2017), dissolved nutrient concentrations were not normally distributed. Therefore, for consistency, we used non-parametric statistics for all comparisons between floodwater and in-estuary concentrations for all analytes measured in all years. This necessitated a reanalysis of the 2017 data but does not invalidate the earlier conclusions (Macias-Tapia et al. 2021). To evaluate the differences between floodwater and baseline estuarine conditions, we used the one-tailed Mann-Whitney rank test. For our study, the alternate hypothesis (H A ) was that the concentrations of samples in the floodwater group were higher than in the group of samples collected in the estuary prior to flooding. If p < 0.05, we accept H A .
To evaluate the variables that controlled the differences in water quality parameters (i.e., dissolved nutrients concentrations and Enterococcus abundance) between floodwater and estuarine pre-flooding water, we (1) analyzed the correlation between the longitude at which the floodwater sample was collected (as a proxy of distance from the mouth of the estuary) and the magnitude of the different water quality parameters and (2) and evaluated the effect of land coverage data from which each sample was collected. For the longitudinal analysis, we used the stats.linregress function, from the SciPy library in Python, to calculate if the linear regression model was significant. When p < 0.05 for the linear regression, we calculated a corresponding R 2 value. For land coverage analysis, we used data from the Virginia Geographic Information Network (VGIN) for the Lafayette River watershed. To group the samples at the land coverage in which they were collected, we used a spatial join function in ArcMap™. We used the stats.mannwhitneyu function in Python to perform a non-parametric signed-rank test to compare the magnitude of the water quality parameters among land coverage groups. When the p-val < 0.05, the differences were considered significant.

Hydrodynamic Model Calculations for the Inundation Volume per Tidal Flooding Event
A more detailed procedure for tidal flooding inundation volumes can be found in Macias-Tapia et al. (2021). Briefly, geospatial inundation depths from a street-level hydrodynamic model, driven in a one-way nested configuration fed by water levels predicted using the Semi-implicit Crossscale Hydroscience Integrated System Model (SCHISM) hydrodynamic model, were used as high-resolution rasters of flood prior to fall king tide events from 2017 to 2021 (Loftis et al. 2018(Loftis et al. , 2019. These inundation depths were used to estimate the tidal flooding inundation volume per event from GIS raster outputs (Danielson et al. 2016

Nutrient Loading Calculation
Nutrient loads during tidal flooding events were calculated using a similar approach to that described in Macias-Tapia et al. (2021). Briefly, we took the differences in dissolved nutrient concentrations between floodwater and baseline samples and multiplied them by the volume of floodwater inundating the perimeter of the Lafayette River. Differences in nutrient concentrations between floodwater and baseline estuarine conditions were only calculated when statistically significant. We compared floodwater nutrient concentrations with baseline estuarine concentrations from (1) transect, (2) time series, (3) CBP, and (4) CBP-year ( Table 2). The minimum nutrient load from flooding was calculated as the difference between the floodwater concentration mean and the estuarine baseline concentration plus one SD. The central input scenario for each analyte was calculated as the difference between the median floodwater and median in-estuary nutrient concentrations. The maximum load was calculated as the difference between the mean floodwater concentration plus one SD and the mean baseline estuarine concentration minus one SD. To account for error in the calculation of floodwater volume, loading was calculated using the average inundation volume calculated each year plus or minus the model's uncertainty for that given year.

Sample Size Validation
Differences in the standard deviations (ΔSD) for each water quality parameter measured were calculated using a randomized sub-sampling method of the entire dataset to determine the approximate number of samples needed to accurately estimate the mean floodwater concentration for each analyte (Krebs 2014). For this calculation, we used the following formula: where "SD n(i) " is the standard deviation (SD) of "n(i)," an accumulated number of random sub-samples of the dataset, and "SD n(i)+10 " is the SD of "n(i)" plus ten random subsamples of the dataset. The calculation was repeated until the maximum number of samples available per analyte in each year was reached. We repeated the randomized subsampling of the dataset process 8 times for each year and water quality parameter and calculated an average ΔSD value from this. The optimum sample number to capture the floodwater concentration variability accurately for each analyte was determined as the sample number at which ΔSD approached zero and remained constant (Krebs 2014).

Interannual Variability in Floodwater Sampling and Inundation Volumes
The extent of flooding varied between years (Table 3). Floodwater inundation volumes from the model predictions for the 2017 to 2021 MTM events fell within a narrow range from 3.54 to 4.55 × 10 7 m 3 , respectively. Model uncertainty changed, based on the meteorological uncertainty of the hydrodynamic model input assumptions, ranging from 3.8 to 4.8%. Using this uncertainty, the standard deviation (± SD) for the inundation model ranged from 0.14 to 0.18 × 10 7 m 3 .
(1) ΔSD = SD n(i) − SD n(i) + 10  Table 4 Dissolved nutrient concentrations (µM) in floodwater and estuarine samples collected prior to flooding over the 5-year study period. Also shown is the number of samples collected (n).
The minimum (Min), maximum (Max), and mean concentrations for each dataset along with the standard deviation (SD) of the mean, the median concentration, and the p-value
From the Chesapeake Bay Program (CBP) data , the mean October/November concentrations for DIN, DON, and PO 4 3− were 16.5 ± 9.0, 19.2 ± 5.7, and 1.2 ± 0.6 µM, respectively (Table 4). However, dissolved nutrient concentrations measured at the CBP water quality station in the Lafayette River in days prior to MTM events in each year (CBP-year) were variable. DIN ranged from 11.8 to 19.8 µM, DON concentrations ranged from 14.4 to 18.7 µM between 2017 and 2021, and PO 4 3− concentrations ranged from 0.96 to 1.8 µM.

Comparisons Between Floodwater and Baseline Dissolved Nutrient Concentrations
Although not normally distributed, dissolved nutrient concentrations in 2017 were closer to a normal distribution. This allowed Macias-Tapia et al. (2021) to calculate a single baseline value and perform a paired analysis between estuarine and floodwater concentrations. However, results for samples collected from 2018 to 2021 were farther from a normal distribution, which required the usage of a different statistical approach. For consistency, we reanalyzed the 2017 data using the same non-parametric unpaired statistical analysis necessary for the analysis of the 2018 to 2021 data. This yielded different results for DIN loading for 2017 than those previously reported in Macias-Tapia et al. (2021) but does not invalidate those results.
Out of the 52 comparisons for dissolved nutrient concentrations and baselines from 2017 to 2021, 58% resulted in floodwater being higher (Table 4). The percentage changes depending on the year, analyte, and baseline considered. In general, floodwater nutrient concentrations were higher than the CBP climatologies and the CBP measurement made prior to the floodwater sampling at  (Table 4). On the other hand, dissolved nutrient concentrations in floodwater were usually lower when compared with cruise transect and time series data from the Lafayette River prior to flooding. For all sampling years, the maximum DIN, DON, and PO 4 3− concentrations in floodwaters were generally much higher than the average and maximum baseline concentrations measured during boat transects, at the time series site, and as part of the CBP monitoring program. The only exceptions were DIN concentrations in 2019 and PO 4 3− concentrations in 2020. In 2017, DIN concentrations in floodwater samples were not statistically higher than any of the baseline datasets when compared using non-parametric statistical tests (Table 4). DON and PO 4 3− concentrations in floodwaters were significantly higher than those measured as part of the CBP water quality monitoring program in October 2017 prior to the MTM event that year (CBP-2017). Floodwater PO 4 3− concentrations were also significantly higher than the October-November CBP climatology. In 2018, there were significant differences between DIN concentrations in floodwaters and those measured at the Lafayette River time series data. Floodwater DON concentrations were significantly higher than both the CBP climatologies and the CBP measurements made in October 2018 (CBP-2018). For PO 4 3− , the floodwater concentrations were significantly higher than the CBP measurements made in the days prior to the king tide in 2018 (CBP-2018). In 2019, DIN concentrations in floodwater samples were significantly higher than those measured at the Lafayette River time series site, and floodwater DON concentrations were significantly higher than CBP climatologies and measurements made just prior to the 2019 MTM sampling (CBP and CBP-2019). In contrast, PO 4 3− concentrations in floodwaters during 2019 were significantly greater than all of the baseline measurements except for those made during the cruise transect. In 2020, floodwater DIN concentrations were significantly greater than those measured on the cruise transect prior to MTM sampling. For DON, floodwater concentrations were significantly higher than concentrations reported for all four baseline datasets. However, PO 4 3− floodwater concentrations were only greater than the CBP-2020 concentrations. In 2021, floodwater concentrations of DIN, DON, and PO 4 3− were significantly higher than all of the baseline datasets available that year.

Differences Between Floodwater and Estuarine Nutrient Concentrations
Excess nutrient concentrations were observed in floodwaters relative to baseline estuarine samples during all years (Table S1). The magnitude of the excess nutrient concentration varied between years, depending on the baseline estuarine data and the model used to make the calculation. The range of excess floodwater nutrient concentrations was estimated using conservative, moderate, and extreme scenarios.
In 2017, DIN differences were not calculated because concentrations in floodwater samples were not statistically higher than any of the concentrations observed in any of the types of baseline estuarine data available (i.e., no excess floodwater DIN). Excess DON concentrations were observed in floodwaters relative to the CBP samples collected from a station in the Lafayette River prior to the MTM sampling event (CBP-2017) and ranged from 8.2 to 34.3 µM, and excess PO 4 3− concentrations in floodwater ranged from below detection (BD) to 11.8 µM using CBP climatology data as the baseline estuarine nutrient concentrations. In 2018, excess DIN concentrations in floodwaters were BD to 42.0 µM using the Lafayette River time series data as the baseline estuarine condition, excess DON concentrations ranged from 32.8 to 76.2 µM, using the CBP-2018 data and the CBP climatology as the baseline estuarine condition, respectively, and excess PO 4 3− concentrations in floodwaters ranged from BD to 4.4 µM using the CBP-2018 data as the estuarine baseline concentration. In 2019, excess DIN concentrations in floodwaters ranged from below detection to 11.4 µM using the Lafayette River time series data as the baseline estuarine concentration, excess DON concentrations in floodwaters ranged from 36.7 to 69.0 µM using the CBP climatologies as the baseline estuarine concentrations, and excess PO 4 3− concentrations in floodwaters ranged from 1.4 to 9.4 µM, using the Lafayette River time series and the CBP climatologies, respectively. In 2020, excess DIN concentrations ranged from 0.1 to 12.4 µM using cruise transect data as the estuarine baseline, excess DON in floodwaters ranged from 6.9 to 96.7 µM, using the boat transect and CBP climatologies as baseline estuarine data, and excess PO 4 3− concentrations in floodwater ranged from 0.2 to 0.9 µM, using the CBP-2020 data as the baseline. In 2021, excess DIN concentrations in floodwaters ranged from below detection to 17.3 µM, using the CBP climatologies and the boat transect datasets, respectively, excess DON concentrations in floodwaters ranged from 16.1 to 112.8 µM, using the boat transect and the CBP climatologies, respectively, and excess PO 4 3− concentrations in floodwater ranged from 0.4 to 4.7 µM, when using the boat transect and CBP climatologies, respectively.

Enterococcus Abundance in Floodwater
Enterococcus abundance in floodwaters ranged from 30 to 24,000 MPN 100 ml −1 over the 5 years ( Fig. 2 and Table S2).
The latter (24,000 MPN) is the upper limit of detection for the method used. For , 2018 of the samples, respectively, were above the upper detection limit for the method. Floodwater samples collected during our multi-year campaign showed that fecal contamination in floodwater samples consistently exceeded the threshold (104 MPN 100 ml −1 ) for recreational waters in Virginia (Fig. 2). In any given year in which samples were collected, > 80% of samples collected were above this threshold.
The abundance of Enterococcus did not correlate with the extent of flooding (Table S3). Enterococcus abundance was highest in 2018 when flooding was modest (0.37 m above MHHW) and intermediate in 2021, when the highest water levels were recorded (0.64 m above MHHW). Enterococcus abundances in 2018 were significantly higher than those observed in any other year. Enterococcus abundances during 2017 and 2021 were statistically similar to each other and statistically higher than abundances on samples collected Table 5 Net input of dissolved nutrients (in × 10 3 kg N or P year −1 ) during annual king tide tidal flooding events. Baseline indicates the preflooding in-estuary condition used to make the calculation (see Table 4). For the different nutrient loading scenarios, we considered a conservative scenario where the mean baseline estuarine concentration (BL) plus 1 standard deviation (1SD) was subtracted from the mean floodwater concentration (FW), a moderate scenario where the median BL concentration (BL med) was subtracted from the median FW (FW med) concentration, and an extreme scenario where the mean BL − 1SD was subtracted from the FW + 1SD. To calculate a total load, we considered the error in our floodwater volume estimate and multiplied the estimated floodwater concentrations by either the mean estimated floodwater volume minus 1SD (Vol low) or the mean floodwater volume plus 1SD. BD indicates that loading for that calculation was below detection during 2019 and 2020. Results from 2020 were only statistically higher than those from 2019.

Sample Size Validation
For most of our sample size validation analyses of floodwater, we found that the change in the standard deviation for dissolved nutrients (ΔSD = SD n(i) − SD n(i)+10 ) decreased and reached or approached zero with the addition of more samples to the analysis in most years (Fig. 3). For DIN, there was more variability in the ΔSD during 2017, 2018, and 2020, but when the sample size was large, the ΔSD reached or approached zero in all years (Fig. 3a). During 2019 and 2021, the values for DIN-ΔSD remained low even at small sample sizes. For DON, the ΔSD for concentrations measured in floodwater also decreased as sample size increased in all years but 2021 (Fig. 3b).
For PO 4 3− , ΔSD values were low and remained constant through the analysis for most years except 2017 (Fig. 3c). Fewer samples were collected to measure Enterococcus abundance (40 per year), and not surprisingly, the ΔSD did not reach zero in any year. However, the ΔSD did decrease as the sample number included in the analyses increased.

Controls on Floodwater Water Quality
There were little to no differences between the longitudinal position in which floodwater samples were collected and the magnitude of the different water quality parameters measured ( Fig. 4A-T). For floodwater DIN concentrations, there was only correlated with longitude in 2017, 2018, and 2021 (R 2 values ~ 0.02, 0.11, and 0.04, respectively) ( Fig. 4A-E). DON concentrations were correlated with longitude only during 2019 and 2020 (R 2 values < 0.1, p < 6.8 × 10 −3 ) (Fig. 4F-J). In keeping with DON and DIN, PO 4 3− concentrations showed a weak correlation with longitude in only 2 of the 5 years (2018 and 2019; R 2 < 0.1, p < 1.8 × 10 −2 ) (Fig. 4K-O). No significant correlation was ever observed between Enterococcus abundance and sample longitude (Fig. 4P-T).
There was no consistency on the statistical differences among land coverage categories for the different sampling years and dissolved inorganic nutrients measured (Fig. 5 and Table S4). For example, floodwater DIN concentrations collected on impervious surfaces were significantly greater than those collected on samples above wetlands from 2017 to 2019, but not for those collected on 2020 or 2021. Furthermore, floodwater DIN concentrations collected on impervious surfaces were significantly greater than those from samples surrounded by turfgrass and trees only in 2017 to 2018, respectively. Similar inconsistencies in the results were found when comparing floodwater concentrations of DON and PO 4 3− among the land coverage categories in which samples were collected.

Discussion
Although there have been studies that point to the relationship between tidal flooding and estuarine water quality (Kiaghadi and Rifai 2019;Smith et al. 2021), this is the first time that the abundance of fecal matter indicators and nutrient loading from tidal flooding events has been quantified in such extensive campaigns, both spatially and temporally. We conducted an inaugural . This work demonstrated the profound effects of tidal flooding on nutrient and fecal material loading to a sub-estuary of the lower Chesapeake Bay, during a perigean spring tide (king tide) flooding event in 2017. Here, we report results describing nutrient loading and fecal material contamination during autumnal king tides spanning 5 years, over which weather conditions and the extent of flooding varied. This work demonstrated that both dissolved nutrients introduced to the estuary during tidal flooding events and the abundance of fecal indicators are beyond local regulations designed to protect water quality and human health ( Fig. 2 and Table 5). Given the increasing frequency and intensity of tidal flooding events in many coastal areas around the world due to sea level rise, these loads need to be considered when setting future restoration targets and measures for protecting human health (Nicholls and Cazenave 2010;Spanger-Siegfried et al. 2014;Neumann et al. 2015;Ezer 2018).

Tidal Flooding Events as Significant Non-point Sources of Dissolved Nutrients
Results from this study demonstrated that dissolved nutrient concentrations were higher in floodwaters than in estuarine waters prior to flooding (Table 4); however, the magnitude of that difference varied between years and was highly dependent on the baseline dataset used to make the comparison. The heterogeneity in the results was likely the result of temporal and spatial variability in estuarine conditions over seasonal, tidal, and diurnal cycles (Morse et al. 2014). Despite this variability in the differences between floodwater and estuarine concentrations, estimates of dissolved nutrient loads during tidal flooding events in all 5 years represented a substantial nutrient load to the estuary (Table 6). We compared our load estimates to the total maximum daily loads (TMDL) for total nitrogen (TN) and total phosphorus (TP) established by US Environmental Protection Agency for TN and TP for the Lafayette River (7.95 × 10 4 and 5.37 × 10 3 kg year −1 , respectively) (EPA 2010). The EPA has set tributary-specific TMDLs for the entire Chesapeake Bay, and these regulations include separate allocations for land-based inputs (e.g., non-point sources such as runoff), point sources including commercial and wastewater discharges, and atmospheric inputs. Given the extensive spatial extension of tidal flooding, we considered our nutrient load estimates from tidal flooding a non-point source of dissolved nutrients and compared them to load allocations for land-based inputs of TN and TP (880 and 58.1 kg year −1 , respectively). We did not measure particulate N in floodwater because of the sampling methods employed but found that the total dissolved N (TDN = DIN + DON) load from floodwater inputs always exceeded the annual land-based load allocation for TN established in the Lafayette River TMDL (Table 6). Our estimates suggest that DIN and DON loads from tidal flooding ranged from 9.5 to 226% of the total load allocation using conservative to extreme scenarios. In the case of phosphorus (P), we also did not measure the particulate fraction, but our results indicate that loads of dissolved inorganic P (DIP) delivered during a single flooding event introduced more P than 100% of the allowable landbased TP load allocations in the Lafayette River TMDL, independent of the scenario considered (Table 6). Our estimates suggest that tidal flooding could introduce 5.1 to over 100% of the total allowable load using both conservative and moderate loading scenarios, but more than 200% when considering the extreme scenario. ) from tidal flooding from the different loading scenarios in Table 5 relative to the Lafayette River TMDL load allocation for land-based inputs (land) and the total load allocation (total). BD indicates that the floodwater concentration was lower than the baseline concentration and so no load could be calculated Even at the more conservative end of the estimations, results from this study show that tidal flooding is a significant source of DIN, DON, and DIP to coastal waterways ( Table 6). The estimates of dissolved nutrients load during tidal flooding events varied up to an order of magnitude by year and baseline used for the calculation. In terms of inorganic nutrients, DIN and DIP, which were large components of the loading we observed, are thought to be more bioavailable to marine microbes than organic N or P (Mousing et al. 2018). Primary production in the Lafayette River, like much of the lower Chesapeake Bay, is thought to be mainly limited by N, and areas in which phytoplankton communities are limited by neither N nor P have decreased during the last decades (Zhang et al. 2021). Dissolved inorganic nutrients added during tidal flooding events could represent another source of nutrients for algal blooms, allowing them to be sustained longer and distributed to other parts of the estuarine system. Prior work has demonstrated that dissolved nutrients concentrations play a major role in algal bloom origination and proliferation in the Lafayette River (Morse et al. 2011(Morse et al. , 2013(Morse et al. , 2014Mulholland et al. 2018) and can promote toxic strains of cyanobacterial blooms more than nontoxic strains (Davis et al. 2010). Furthermore, our results indicate that the percentage of P loads during a single tidal flooding event, relative to the TMDLs, is higher than those for N (Table 6). This stoichiometric imbalance affects phytoplankton community growth rates, species composition, and nutritional quality (Arrigo 2005;Burson et al. 2016).
Previous studies have shown that some bloom-forming algal species are mixotrophs (e.g., Margalefidinium polykrikoides) and could also benefit from DON inputs (Mulholland et al. 2018). We did not measure dissolved organic P (DOP), but previous studies have pointed to higher DOP uptake rates by phytoplankton when DON is introduced (Fitzsimons et al. 2020). Additionally, there are other potential contaminants that could be introduced to estuaries as a result of tidal flooding that play a role in the water quality of coastal systems. These include other nutrients, like Mg, Na, and S (Hofmann et al. 2021), heavy metals (Rakib et al. 2022;Wu et al. 2020;Price et al. 2019), microplastics (Dodson et al. 2020), and contaminants of emerging increasing concern, including pharmaceuticals and personal care products (Kroon et al. 2020). Future studies should assess how these inputs may be changing as a result of increases in tidal flooding and their influence on human health. Failing to account for increased non-point source loading of nutrients and contaminants due to increases in tidal flooding could result in further deterioration of the systems we are trying to restore.

Controls on the Magnitude of Nutrient Loadings from Tidal Flooding
We initially hypothesized that nutrient loading from tidal flooding would correlate with the extent of inundation during flood tides. However, results from our multi-year analyses show that there were no significant correlations between the floodwater volume and the nutrient load delivered during the five annual events included in this study (Fig. 6). Further analysis showed that there was no correlation between nutrient or Enterococcus concentrations with the distance from the mouth of the estuary at which the floodwater sample was collected (Fig. 4) or with the land coverage characterizing the sampling point (Fig. 5). In terms of land coverage, the Lafayette River watershed is highly residential (Fig. S1). This makes the distribution of these properties very homogeneous among floodwater samples and probably hides the effect of land use categories on dissolved nutrient concentrations. Multiple studies have shown that there is a relationship between land use and nutrient runoff (Tu 2011), with agricultural and urban areas delivering the largest nutrient loads and exerting the most water quality impacts (e.g., Cheng et al. 2022). However, these studies integrated land use over a much larger area than our small sub-watershed and that floodwater nutrient loads may need to be similarly integrated over larger spatial scales. Future studies should consider sampling multiple watersheds wherein the overall land use and land cover are more distinct (e.g., agriculture vs industrial). Weather conditions and landscape prehistory before and during flooding can also modify the magnitude of nutrient inputs. For example, precipitation and runoff prior to flooding can remove materials from the landscape (Selbig 2016), which will attenuate nutrient inputs during subsequent flooding. Sampling for this study was focused only around perigean-spring tides occurring during the fall season. This resulted in homogeneous weather conditions during the 5 years of our sampling campaign. Future studies are needed to quantitatively assess the temporal variability of nutrient inputs during tidal flooding with respect to the immediate prehistory of the landscape.
Apart from the critically high nutrient inputs associated with tidal flooding events, relative to establish restoration limits, results from our study indicate that there are hotspots where concentrations of dissolved nutrients were orders of magnitude higher than estuarine concentrations measured prior to flooding (Table 3). In all years, outlier concentrations were found in samples collected at or near public parks. Several studies indicate that in green areas located in or near cities, the presence of pets can elevate N and P concentrations (John et al. 2020;Paradeis et al. 2013;De Frenne et al. 2022). Our study suggests that interactions between floodwaters and green areas designated for pets were hotspots for nutrient loading. Other hotspots that might have contributed disproportionately to nutrient loading from tidal flooding included sites near a stormwater retention pond and residential areas with a large areal proportion of impervious surfaces. A recent study found that the design and age of stormwater retention ponds can modify their ability to retain nitrogen (Gold et al. 2017). Other studies suggest that runoff from residential and urban areas contributes substantial loads of DIN, DON, and DIP to connected waterways (Toor et al. 2017;Lian et al. 2019;Jani et al. 2020). Our work indicates that identifying and targeting hotspots of nutrient inputs from flooding events will be critical for prioritizing management actions aimed at minimizing nutrient loading from tidal flooding.

Enterococcus Abundance in Floodwater
Floodwater samples collected during our multi-year campaign showed consistently elevated Enterococcus abundances in floodwater samples, with values well above the threshold for safe recreational use of water (104 MPN 100 ml −1 ) (Fig. S2). Ezer (2018) found an increase in the number of hours per year in which the water level in the city of Norfolk is 0.53 m above the mean higher high water (MHHW). This is worrisome as humans and animals are in intentional and unintentional contact with floodwaters putting their health at risk (Iversen et al. 2004). Therefore, identifying the sources of fecal contamination and hotspots where these contaminants accumulate is critical for the health and resilience of urban areas that experience recurrent tidal flooding. Hart et al. (2020) pointed to the interaction between floodwater and stormwater infrastructure (Hart et al. 2020), while data-driven models built with high frequency sampling indicated that the most important variables to predict fecal indicator bacteria were tide, solar irradiation, water temperature, significant wave height, and offshore wind speed (Searcy and Boehm 2021).
Our multi-year analysis indicated that there was no statistical correlation between the abundance of Enterococcus and the extent of flooding (Table S3). Other factors, like land coverage and use, can influence the abundance of fecal indicators in the water. Such factors include runoff from high-and low-density residential areas, such as those found in the Lafayette River watershed (Tiefenthaler et al. 2011),  Table S5 and water in contact with drain systems (Haile et al. 1999;Aryal et al. 2021). Future studies should include a higher number of samples that represent all of the affected areas and a more conclusive land use characterization system. Furthermore, samples representing pre-flooding conditions should be considered, similar to the analysis done with the dissolved nutrients in this study.

Variability Analysis: Successful Floodwater Characterization with a Citizen Science Approach
Capturing the variability of dissolved nutrient concentrations in floodwater is important for quantifying nutrient loads at a watershed level. This becomes more relevant when considering the highly variable biogeochemistry of estuarine systems, like the Lafayette River (Morse et al. 2014) that are affected by multiple biotic (Mulholland et al. 2009(Mulholland et al. , 2018 and abiotic factors (Morse et al. 2011(Morse et al. , 2013. Although further analysis suggests that the number of samples collected and sampling events performed were not enough to capture the variability of the system as a whole (Fig. 7), our results demonstrate that citizen science sampling at a watershed scale can provide enough data to accurately quantify nutrient loading from individual tidal flooding events by evaluating the variability in SD (i.e., ΔSD) when increasing the number of random floodwater samples collected. Although citizen science projects have proven to be cost-effective ways to gather data (Tulloch et al. 2013;Ferri et al. 2020), organizing and training volunteers to collect floodwater sampling takes monetary and human resources (Alfonso et al. 2022). This study can advise future work by providing an estimate of how many samples are needed to produce statistically robust estimates of nutrient loading due to tidal flooding. Our results confirmed that a large sample size (> 100 samples) is needed to fully capture the variability in all of the targeted analytes at the scale of the Lafayette River watershed (Fig. 3) and robustly estimate nutrient loading from tidal flooding (Tables 3, 4, and 5). Future studies are needed to constrain the variability in nutrient loads due to other factors that might contribute to variability in excess nutrients in floodwaters (e.g., variability in landscape condition just prior to flooding, time since last rainfall, seasonality) so that we can better estimate nutrient loading from tidal flooding and reduce these loads.

Conclusions
Sea level is rising, leading to an increase in tidal flooding in low-lying coastal communities like our study site (Ezer 2018). Despite widespread recognition of this trend and the economical effects it has, little is known about nutrient loads delivered during tidal flooding and its impact on the water quality. In this study, we quantified nutrient loading to a lower Chesapeake Bay sub-tributary, during king tides from 2017 to 2021. We leveraged a citizen science campaign to collect a sufficient number of samples at high tide. We found the following: • DIN, DON, and PO 4 3− concentrations were higher in floodwaters than in estuarine waters, but the magnitude of this difference was highly dependent on the baseline estuarine dataset used as a reference. • Loads of DIN, DON, and PO 4 3− from singular flooding events can contribute more than 100% of the annual tributary-specific total and land-based federal allocations (i.e., total maximum daily loads (TMDL)) for total nitrogen (TN) and total phosphorus (TP). It is important to notice that the magnitude of the estimated load varied up to an order of magnitude depending on the year and pre-flooding conditions considered. • Enterococcus abundances in floodwater usually exceeded the threshold for recreational water use in Virginia (104 MPN 100 ml −1 ) during the 5 years of sampling. • Results from our temporally and spatially extensive sampling campaigns showed that there were no significant correlations between Enterococcus abundance or dissolved N and P concentrations in floodwater and (1) the inundation volume in a single event, (2) the distance from the mouth of the estuary at which the floodwater sample was collected, or (3) the land coverage characterizing the sampling point. • Our analysis of the variability of dissolved nutrient concentrations in floodwater indicates that > 100 samples are needed to produce statistically robust estimates of nutrient loading due to tidal flooding. Using our citizen science approach, we met that goal.
Results from this study indicate that tidal flooding contributes a substantial nutrient load to the aquatic systems we are trying to restore. This load needs to be included in management and policy decisions aimed at managing nonpoint loads of nutrients. Additionally, our results indicate that floodwater events represent a risk for human health due to high abundance of fecal indicators. Measures should be taken to control this or to inform people of the risk.