Biogeochemical control points of connectivity between a tidal creek and its floodplain

As global climates shift, coastal systems experience changes that alter function within the tidal zone. However, it remains uncertain how changes in tidal extent and magnitude will alter coastal biogeochemical cycling. We present high‐frequency data collected in situ along two transects across a tidal creek and floodplain to capture how vertical and lateral connectivity changes act as biogeochemical control points. Using wavelet coherence analysis, we identified short periods of vertical tide–floodplain connectivity during spring tides followed by extended periods of floodplain isolation from tides. During a drier period, neap tides were linked to lateral creek–floodplain connectivity, detected by discharge of saline, oxygen‐depleted groundwater to the creek on ebb tides. We estimate that up to half of creek water throughout the tidal cycle was associated with floodplain drainage. From point measurements of carbon, nitrogen, and greenhouse gases along the creek–floodplain transect, we infer that lateral connectivity functions as a biogeochemically significant control point.

Coastal ecosystems influence global biogeochemical cycles as regulators of carbon, nutrient, and pollutant exchanges among land, ocean, and atmosphere. For example, coastal carbon sequestration rates are an order of magnitude higher than inland systems (Mcleod et al. 2011;Duarte et al. 2013), where estuaries and wetlands account for a disproportionately high percentage of carbon sequestered relative to coastal margins (Najjar et al. 2018). In addition, coastal systems provide a wide range of other ecosystem services, including erosion control, fisheries habitat, and flood mitigation for growing coastal populations vulnerable to rising seas (Barbier et al. 2011;Kulp and Strauss 2019).
Within coastal regions, the intertidal zone functions as a biogeochemical control point (i.e., a point in time and/or space of disproportionate biogeochemical significance, per Bernhardt et al. 2017) that dictates mixing of marine and freshwater inputs and lateral exchange along the terrestrialaquatic continuum (Nezlin et al. 2009;Regier and Jaffé 2016;Bogard et al. 2020). As global climate shifts, factors including sea level rise, altered precipitation regimes, and increased storm activity are changing tidal extent and coastal geomorphology (Ensign and Noe 2018). Concurrently, higher anthropogenic nutrient loads delivered to coastal systems drive increased eutrophication and degraded water quality (Doney 2010;Seitzinger and Phillips 2017).
Together, these changes form a compounding disturbance to coastal physical, chemical, and biological processes (Ward et al. 2020;Yabusaki et al. 2020). A strong, mechanistic understanding of the patterns and drivers of these processes is required to accurately predict how such shifts in biogeochemistry will affect coastal ecosystems, the services they provide, and the populations that depend on them. This is particularly true for small tidal watersheds, which dominate the coastal land-water interface, and are disproportionately vulnerable to sea-level rise (Tagestad et al. 2021). However, because the function of the intertidal zone is highly complex and dynamic, predictive modeling of these processes remains difficult, leaving potential future impacts largely unconstrained (Ward et al. 2020).
In this study, we explored connectivity and resulting biogeochemical control points at the terrestrial-aquatic interface in the tidal zone. We collected high-resolution hydrology and water quality data in a tidal creek and its floodplain to identify periods of connectivity. Next, we used wavelet coherency analysis, a tool for visualizing time-and frequency-varying relationships (Carey et al. 2013;Alexander et al. 2020) to identify periods of connectivity between tidal flooding and groundwater, and a simple linear regression model to quantify groundwater discharge back into the tidal creek. We used these periods to establish a conceptual model that relates hydrology and water quality measurements to control points of coastal ecosystem function and infer biogeochemical importance.

Site description
We established a network of aquatic sensors at Beaver Creek (BC), a first-order stream in Washington State. BC drains a catchment of 382 ha, and is characterized by forested hillslopes that bound lower elevation floodplains, which account for 13% of the catchment area (Yabusaki et al. 2020). BC is a tributary to the Johns River near its confluence with Grays Harbor, which flows into the Pacific Ocean (Fig. 1). In 2014, Washington Fish and Wildlife removed a barrier to saltwater intrusion (WDFW 2019), reexposing the lower reaches of BC and associated floodplains to semidiurnal tides.
Data used in this study were collected from 2018 to 2019 at sites located mid-channel in the creek (C down and C up ) and in 1 m deep, 5.1 cm diameter shallow groundwater wells located at the confluence of the hillslope and floodplain (W down and W up ) along two transects (downstream and upstream) (Fig. 1). Additional data for floodplain groundwater chemistry were collected at two floodplain sites between C up and W up (F up1 and F up2 ) (Fig. 1). C down experiences stronger tidal influence than C up (mean salinity is 2.5 PSU [ 35%] higher) although both sites experienced salinities > 28 PSU during the drier summer months. Unpublished data collected in wells located in the upstream floodplain indicate strong reducing conditions, high oxygen demand, and persistent anoxia in floodplain groundwater.
We discuss two forms of connectivity between surface waters in the creek and groundwater in the floodplain. Vertical connectivity refers to tidal flooding when creek water infiltrates vertically into floodplain soils. Lateral connectivity refers to groundwater in the floodplain moving laterally into the creek (Fig. 1).

Data collection, correction, and analysis
C down and C up were equipped with EXO2 water quality sondes (YSI) measuring temperature, pressure, conductivity, dissolved oxygen (DO), pH, and turbidity. W down and W up were equipped with HOBO sensors (Onset Corporation) measuring temperature, pressure, and conductivity. Sensors were calibrated per manufacturer protocols during maintenance visits and measured at 30-min intervals. Differences between pre-and post-cleaning readings were used to correct offsets associated with fouling/drift. Hydraulic head, which equalizes water levels between sites to a common vertical datum, and salinity were calculated from sensor data following methods in Yabusaki et al. (2020).
Rainfall was measured using an Onset rain gauge near F up1 . We selected two time periods to analyze, which we refer to as Dry (12 July 2018 to 01 October 2018) and Wet (01 February 2019 to 03 May 2019). These time periods were selected for two reasons. First, they have high quality data at all sites over two consecutive spring-neap cycles, allowing us to investigate the impact of lunar-scale tidal cycles on creek-floodplain connectivity. Second, they represent distinct periods of seasonal hydrology, observed both by differences in rainfall, and differences in surface and groundwater salinity.
Point samples were collected during sensor maintenance visits for additional biogeochemical context. Samples for dissolved organic carbon and total dissolved nitrogen were filtered (0.2 μm) then analyzed using a Shimadzu TOC-L analyzer equipped with a total nitrogen module. Samples for pCO 2 and pCH 4 were collected following Ward et al. (2019), and analyzed by direct injection on a Picarro G2508 Cavity Ring-Down Spectrometer. Gas samples were collected into pre-evacuated gas-tight glass vials and analyzed within 1 week of collection. Additional details on point sampling are presented in Supplemental.

Statistics
Statistical analyses were conducted in R version 4.0.2 (R Core Team, 2020). We used Wilcoxon tests to compare means and lm functions for regression models. To estimate lateral inputs to the creek during periods with decreased DO concentrations (DO sags) at C up (Fig. S1), we constructed linear regression models for salinity/DO between C up and C down Grays Harbor is provided for geographical context in the upper right. Elevation profiles of upstream and downstream transects show the location of each site laterally and vertically relative to the tidal creek. Circles represent sites with sensors, and triangles denote sites where point samples were collected, but no sensors were deployed. The shaded blue area represents the maximum head values for C down and C up during the time periods presented in Fig. 2. Both creek sensor site abbreviations start with "C" while both well sensor site abbreviations start with "W." The downstream transect site abbreviations all contain "down" subscripts, while the upstream transect site abbreviations all contain "up" subscripts.
during tides without DO sags ( Fig. S2) to establish baseline relationships. We then used these relationships and salinity/ DO time series at C down (which did not experience DO sags) to predict values for concurrent tides with DO sags at C up (Fig. S2). Finally, we estimated the percentage of groundwater inputs as the average of residuals (observedpredicted) for each tide at C up using salinity and DO as independent predictors (Fig. S3). See Supplemental for a full description of the approach and assumptions made.

Wavelet coherence analysis
Wavelet coherence analysis applies a continuous wavelet transform to bivariate time series to analyze the coherence (localized correlation coefficients) between two variables in the time-frequency domain (Roesch and Schmidbauer 2018), and has been used to interpret inter-variable relationships in environmental datasets (Carey et al. 2013;Alexander et al. 2020). This approach isolates coherence between two variables at different frequencies (e.g., subtidal, tidal, daily, weekly, and spring-neap), which is useful for noisy signals like tidal hydrology and water quality time series. In our coastal system, these frequencies are largely tied to cycles (tidal, diurnal, neap-tide), so we label frequencies as cycles. We used the WaveletComp R package (Roesch and Schmidbauer 2018) to conduct wavelet coherence analysis based on a mother Morlet wavelet for head and salinity at the two well sites (W up and W down ) to examine when hydrology (represented by head) and water chemistry (represented by salinity) are locally correlated. We interpret coherence between these parameters as "vertical connectivity" where tidal flooding is increasing salinity at the frequency of a given cycle.

Comparing wet and dry periods
During the Dry period, rainfall occurred as isolated events, averaging less than 1 mm d −1 (Fig. S4A), compared to more consistent rainfall throughout the Wet period, averaging almost 4 mm d −1 (Fig. S4B). Consequently, hydraulic head was lower across all sites during the Dry period (Fig. 2), allowing stronger tidal exchange between the creek and unsaturated floodplain soils.
During both Wet and Dry periods, the floodplain hydraulic head was influenced by tides in connection to spring-neap tidal cycles (Fig. 2). However, when comparing Dry and Wet time series, there is a clear difference in well salinity patterns, which increased during Dry spring tides (i.e., mid-July, mid-August, and mid-September), but did not increase during Wet spring tides. Similarly, salinity and head both decreased during Dry neap tides, but neither decreased notably throughout Wet neap tides (Fig. 2). During the Dry period, salinity at W up is higher than W down during spring tides and lower than W down during neap tides, but during the Wet period, salinity at W down is consistently 2 units higher than W up .
Wavelet coherence: Vertical connectivity of tides to the floodplain We applied wavelet coherence analysis to explore the relationship between head and salinity at the well sites during Dry and Wet periods (Fig. 3). This relationship indicates the influence of tidal flooding on groundwater in the floodplain, where variation in well head at tidal timescales represents physical inundation, but not necessarily groundwater recharge, while

Regier et al.
Tidal creek biogeochemistry tidal variation in well salinity indicates groundwater salinization (i.e., tidal groundwater recharge). Thus, we interpret stronger coherence between head and salinity at tidal timescales as indicative of stronger vertical connectivity between surface flooding and floodplain groundwater. During the Dry period, strong coherence between head and salinity is present at diurnal and tidal timescales for short periods at both well sites, indicating reoccurring, ephemeral vertical connectivity between creek and floodplain (Fig. 3). Dry period connectivity occurred during spring tides, consistent with Fig. 2, and appeared to be stronger at W up than W down , with no clear coherence at diurnal or tidal timescales during the September spring tide at W down (Fig. 3). The difference in average power between W up and W down for the Dry period indicates stronger coherence at W up for weekly and diurnal cycles (Fig. S5). During the Wet period, we observed stronger coherence at W up again for weekly and diurnal cycles, but also at tidal cycles (Fig. S5), as visually observed by stronger spring tide vertical connectivity at W up relative to W down (Fig. 2).
Water quality excursions: Lateral connectivity of the floodplain to the creek During the Dry period, we observed two phases where DO sagged to hypoxia then recovered to equilibrium during several consecutive tidal cycles at C up , but not C down (Fig. S1).
Consistently, when DO values are lower than expected, salinity is higher than expected (Fig. S2). Because sags occur on the ebb tide, the source of these water quality excursions is not estuarine. Furthermore, because upstream surface inputs would be oxygen-saturated, less saline, and minimal during the Dry period, we attribute the source of water quality excursions to groundwater moving laterally from the floodplain back into the creek. The chemistry and timing of suggested lateral inputs are consistent with similar reports from different coastal systems, including other tidal creeks, marshes, and mangroves (Mattone and Sheaves 2017;Nelson et al. 2017;Trifunovic et al. 2020).
Based on this explanation, we used regression models to estimate the percentage of creek water during DO sags attributable to groundwater. We used the values for average DO decrease and salinity increase for each tide analyzed (Fig. 4A), averaging −2.24 mg L −1 and 2.30 PSU, respectively, to estimate what percentage of each tide is attributed to groundwater inputs (Fig. 4B). Based on our assumptions (see Supplemental), decreases in DO suggest that groundwater inputs account for approximately 47% of tidal volume, while increases in salinity suggest approximately 27% (Fig. 4B).

Spatial floodplain groundwater chemistry
Point samples collected at C up , F up1 , F up2 , and W up (Fig. 1 creek and floodplain (Fig. 4C) and provide chemical context for sensor observations. Consistent with our proposed explanation of lateral connectivity, DO is significantly lower at all floodplain sites compared to C up . Likewise, we observed significantly higher values for dissolved organic carbon, total dissolved nitrogen, and pCO 2 , and pCH 4 in the floodplain relative to the creek.

Drivers of vertical tide-floodplain connectivity
We identified three potential drivers of observed vertical connectivity: spring-neap cycles, rainfall, and floodplain geomorphology. First, coherence only occurred at peak salinity during spring tides (Fig. 3), indicating that the lunar calendar is a primary driver. As spring tides yield minimal change in well salinity during the Wet period compared to the Dry period (Fig. 2C), our data suggest the seasonality of rainfall as a second driver, where reduced freshwater inputs would lower the floodplain water table, creating a deeper vadose zone, facilitating groundwater infiltration during tidal inundation. We tested this hypothesis by examining the relationships of freshwater and tidal inputs (using rainfall and well salinity as proxies, respectively) to vadose zone depth (using well head as a proxy) and found a significant correlation to freshwater, but not tidal inputs, supporting our hypothesis (Fig. S6).
Third, the unexpected finding that W up exhibits stronger coherence between head and salinity during spring tides compared to W down (Fig. 3) suggests an additional geomorphic control of vertical connectivity. The average gradient from floodplain to creek at W down is less than half as steep as W up (0.011 m rise per m run and 0.027 m rise per m run, respectively). A shallower gradient and greater distances between W down and C down relative to the upstream transect fits slower drainage, potentially limiting vadose zone depth at W down relative to W up . It is also possible that other geomorphic features (soil density, macropores, etc.) contribute to differences in floodplain drainage.
Together, high-resolution time series of head and salinity, and the relationship between them derived from wavelet coherence, help untangle the roles of seasonal rainfall, tidal  Fig. 1) for salinity, DO, dissolved organic carbon (DOC), total dissolved nitrogen (TDN), and partial pressures of carbon dioxide (pCO 2 ) and methane (pCH 4 ). Asterisks indicate a significant difference in mean between a given floodplain site and the creek based on Wilcoxon tests (p < 0.05).
cycles, and geomorphology in controlling vertical connectivity between the tidal creek and its floodplain. We note that wavelet coherence is particularly useful for this type of analysis, where the periodic coherence of hydraulic head and salinity is made clear, even during the Wet period (Fig. 3), while no covariation is observable in the raw time series (Fig. 2).
Our results underscore wavelet analysis as an easy-to-use yet powerful tool for understanding the importance of distinct temporal cycles in dynamic ecosystems (Nezlin et al. 2009).
Lateral floodplain-creek connectivity as a coastal control point DO sags coinciding with the end of neap tides, when creek water levels are at their lowest, and only observed at the upstream creek site, suggest an ecosystem control point (Bernhardt et al. 2017) isolated both in time and space. Following vertical connectivity during spring tides (Fig. 3), as inputs infiltrate into the floodplain, biotic and abiotic (redox geochemistry) processes remove oxygen (Yabusaki et al. 2017). During the neap tide that follows, groundwater inputs return saline, oxygen-poor water back into the creek (Figs. S1-S3; Fig. 4). Because there are multiple potential sources of salinity and oxygen along this transect, either DO or salinity may act as a more conservative tracer for groundwater inputs in our system. Thus, we estimate lateral inputs from the floodplain account for approximately a quarter (based on salinity) to a half (based on DO) of the water in the creek during Dry period neap tides at W up (Fig. 4B). Our estimates are more conservative than estimates based on 222 Rn as a tracer in a different tidal creek, which suggested that 76% of creek surface water was groundwaterderived, resulting in high CO 2 emissions (Atkins et al. 2013).
As tidal inputs move through the floodplain toward the creek, they come in contact with (and potentially become enriched in) carbon, nutrients, and greenhouse gases (i.e., CO 2 and CH 4 ). Although most point samples do not overlap temporally with DO sags, spatial gradients across the floodplain suggest the chemistry of lateral inputs sourced from floodplain groundwater. We therefore suggest the lateral connectivity events identified by sensor data represent "export" control points (Bernhardt et al. 2017), where hydraulic gradients drive short, rapid, and biogeochemically significant pulses, as is observed in other tidal systems (Regier and Jaffé 2016;Bogard et al. 2020) and soil incubations performed at Beaver Creek (Sengupta et al. 2021). Carbon and nutrients mobilized from groundwater into the creek may represent outwelling events that support productivity of downstream receiving waters, or, alternatively instigate eutrophication events downstream (Seitzinger and Phillips 2017;Bogard et al. 2020;Tagestad et al. 2021), with potential impacts to coastal oceans (Doney 2010). GHGs mobilized via lateral inputs into the creek, as observed in other tidal systems (Sadat-Noori et al. 2016;Trifunovic et al. 2020), may reduce the efficacy of the floodplain to sequester carbon, while creating an additional source of GHGs. Because equivalent decreases/increases in DO/salinity on ebb tides were not observed at C down (Fig. S1), supported by Dry period transit times between C up and C down estimated from velocity measurements at > 36 h, our data indicate that the isolated, biogeochemically complex events associated with periodic lateral connectivity between the floodplain and creek represent the potential for increased biogeochemical activity, but appear to be highly localized both in time and space.

Implications for coastal biogeochemistry
Taken together, high-frequency data suggest spring tidal cycles drive vertical connectivity of tide to floodplain (Figs. 2,  3), followed by neap tide lateral connectivity (Figs. S1-S3), which pumps water from floodplain sediments rich in carbon, nutrients, and GHGs (Fig. 4) into the creek, and potentially to the estuary, as summarized in a simple conceptual model (Fig. 5). Small tidal creeks (1 st to 4 th order) like BC represent 66% of total tidal stream length along the conterminous US coastline (Tagestad et al. 2021). If the mechanisms in Fig. 5 are present in other tidal creeks, the resulting export of carbon, nutrients, and GHGs may be a globally significant control point for coastal biogeochemistry. As small-order tidal streams are predicted to be more impacted by sea level rise than larger tidal rivers (Tagestad et al. 2021), this control point may shift as coastal systems respond.
Sea level rise projections predict tidal regimes will increase the elevation of high tides, and shift tidal extent inland (Ensign and Noe 2018), increasing floodplain area and pumping saltwater into previously freshwater areas. Simultaneously, climate models predict extended drought (Trenberth 2011), where decreased freshwater inputs (analogous to our Dry period) may create a deeper vadose zone. These predicted changes favor stronger vertical connectivity between tidal channels and their floodplains, although, as identified in our data for a tidal creek and its floodplain. Following vertical connectivity events, oxygen (O 2 ) levels are drawn down by biotic and abiotic processes. As groundwater moves laterally back toward the creek, it interacts with, and may mobilize, considerable pools of carbon (C), nitrogen (N), and greenhouse gases (GHG). discussed above, this appears to be contingent on site-specific geomorphology.
The implications of rising seas and shifting precipitation regimes on lateral connectivity are less clear. If the vertical extent of the tidal zone moves up, this may mean neap ebb tides identified as biogeochemical control points will not drain floodplains as efficiently. However, the lower average water table during the Dry period suggests reduced freshwater inputs resulting from longer periods of drought may offset tidal increases to the water table, maintaining neap-tide groundwater discharge. Under this scenario, we would expect enhanced flushing of floodplains relative to current conditions. We note here that our simple model is based on a macrotidal system with low hydraulic conductivity soils. For mesotidal or microtidal systems, shallower hydraulic gradients between creek and floodplain could reduce flushing rates, potentially reducing efficacy of both vertical and lateral connectivity mechanisms summarized in Fig. 5, whereas more hydraulically conductive soils could increase connectivity.
As coastal ecosystems and the services they provide face a shifting climate, there is growing need to constrain biogeochemical processes at the coastal interface to forecast their future (Ward et al. 2020). However, sufficient data to represent the spatiotemporal variability within and between coastal ecosystem types, and the resulting "big data" resources, require novel approaches to manage and analyze (Durden et al. 2017). Simple, statistically robust tools like wavelet analysis or cumulative sums (Nezlin et al. 2009;Regier et al. 2019;Alexander et al. 2020) can help simplify and automate initial interpretation of large coastal datasets. Combined with growing repositories of publicly available, high-resolution data, such tools can help test site-specific mechanisms (Fig. 5) at other sites, and potentially scale up to regional or continental contexts.