Floodplain Inundation and Salinization From a Recently Restored First‐Order Tidal Stream

The systematic response of coastal ecosystems to inundation and salinity exposure is fundamental to their ecology and biogeochemical function. Here we observe and model freshwater‐seawater interactions in a first‐order stream—floodplain system where tidal access was recently restored. Subsurface flow and transport modeling were used to quantify and better understand the interplay of processes, properties, and conditions that control water level and salinity in the floodplain to the tidal stream. Water levels in the stream were highly correlated with tidal forcing, which resulted in episodic inundation of the floodplain at quasi‐monthly frequency. The tidal stream is the only source of salinity to the floodplain, yet shallow groundwater salinity was considerably higher than average stream salinity. The low‐permeability clay floodplain soils limit lateral groundwater flow and transport, resulting in floodplain groundwater and salinity dynamics driven almost exclusively by infiltration during inundation events. As inundation occurs during high tide, estuarine waters reach the floodplain with minor attenuation in salinity from the stream's freshwater discharge. Infiltration and salinity exposure are topography controlled and regulated by ponding depth and duration, seasonal ground saturation, and depth to water table. The model suggests that floodplain salinity is currently in an early stage of transition from pre‐restoration freshwater conditions and will not reach equilibrium for ~20 years. These findings have broad relevance for understanding how and over what time scales coastal ecosystems will respond to increasing seawater exposure from sea level rise, ocean‐originating storms, and changes in natural and man‐made barriers.


Introduction
Coastal ecosystem dynamics are intricately tied to atmospheric, terrestrial, and marine hydrologic forcings, such as rainfall, freshwater discharge, and tidal cycles, respectively. Thus, these ecosystems represent a nexus between marine, freshwater, and terrestrial processes. Coastal ecosystems, such as marshes, mangroves, tidal rivers, and mudflats, are diverse in their physical structure, hydrologic regimes, and ecological composition (Crain et al., 2004;Steel & Milliken, 2013). Such ecosystems are ecologically and biogeochemically important and serve as important reservoirs of organic matter (Duarte et al., 2013;Spivak et al., 2019). The interplay between the various hydrological forcings in these ecosystems are fundamentally important for ecosystem processes and function (e.g., biogeochemical cycles, vegetation zonation; Atkins et al., 2013;Krause et al., 2017;Moffett et al., 2012). Coastal biomes, such as tidal marshes, display dynamic groundwater, surface water, and porewater interactions dependent on topography and hydraulic properties (Moffett et al., 2012). Coastal aquatic systems, such as tidal rivers, are likewise influenced by the interplay between tidal and fluvial energy and geomorphology (Ensign et al., 2013). Thus, the dominant hydrologic forcings present in coastal ecosystems impact not only their physical structure but also their ecosystem function (Krause et al., 2017). Given that coastal ecosystems are also experiencing change at increasing frequency and intensity, such as from sea level rise (SLR) and ocean-originating storms (Craft et al., 2009;Ensign & Noe, 2018;Hopkinson et al., 2008), understanding the hydrologic dynamics and interplay between freshwater and tidal forcings on coastal terrestrial-aquatic connectivity is critical for predicting their vulnerability.
Modeling and field-based efforts have elucidated drivers of groundwater hydrological flow in complex, coupled terrestrial-aquatic systems such as marsh-river systems and freshwater river-floodplain systems. Drivers of tidally-influenced marsh-river systems include seasonal precipitation and evapotranspiration, upland groundwater dynamics/exchange, and seawater intrusion (Peterson et al., 2019;Wilson & Morris, 2012;Xiao et al., 2017). This complex hydrology impacts the terrestrial ecosystem ecology (e.g., marsh plant zonation), nutrient supply, and redox regimes in coastal systems (Marani et al., 2006;Wilson & Gardner, 2006;Xiao et al., 2017). However, the dominant mechanism(s) behind such surface water-groundwater interaction is system-dependent, with tidal forcings, topography, and hydraulic conductivity emerging as primary controllers of terrestrial-aquatic connectivity in tidal marshes (Moffett et al., 2012;Xin et al., 2011Xin et al., , 2013. Freshwater river-floodplain connectivity is dominated by the two-way exchange of water and nutrients with complex hydrologic conditions (Covino, 2017) and resultant controls on biogeochemical dynamics (Harvey et al., 2019;Tockner et al., 1999). Except in watersheds with extensive lake coverage, anthropogenic irrigation and/or strong evaporation due to arid climates, river connectivity to freshwater floodplains/wetlands is dominated by riverine forcings (Yamazaki et al., 2011). Terrestrial-aquatic connectivity in these coupled systems are thus largely driven by flood pulse events (Junk et al., 1989), increasing complexity of surface-groundwater flow and interactions (Khublaryan & Zyryanov, 2006). Such emergent hydrologic drivers presiding over marsh-river systems and freshwater river-floodplain systems are likely present in tidal river-floodplain systems, as tidal river reaches experience nonstationary tides resultant of channel morphology and river flow (Buschman et al., 2009;Matte et al., 2013); however, the hydrologic connectivity between tidal river-floodplain systems is comparatively understudied. For example, the extent and duration of flooding in tidal river floodplains may be daily, seasonal, or episodic and depends largely on elevation and geomorphology of the river system (Ensign et al., 2014).
Tidal inundation is a defining hydrologic component of coastal ecosystems, where seawater floods the landscape on tidal cycle intervals. Coastal floodplains subject to tidal inundation are common occurrences, with relative SLR differentially impacting coastal regions depending on local geological and anthropogenic forcing (e.g., plate subduction, uplift, isostatic rebound, sediment compaction, oil and gas extraction) (Hopkinson et al., 2008). These floodplains and their connectivity to tidal river systems are likely to be increasingly impacted by SLR (Ensign & Noe, 2018); in fact, coastal inundation increases when tidal range is considered. Since around 0.4% of the contiguous United States (~32,000 km 2 ) lies less than 1 m above the mean high water line (MHWL) (Strauss et al., 2012), SLR will increase the extent and duration of coastal inundation, which increases terrestrial exposure to seawater and connectivity between tidal rivers and the surrounding landscape.
Salinity incursion into freshwater wetland reaches of tidal rivers is seasonally variable (Anderson & Lockaby, 2012), and tidal influences on coastal terrestrial systems are linked to elevation and geomorphology (Dalrymple & Choi, 2007;Ensign et al., 2013;Hoitink & Jay, 2016;Langley et al., 2009). Salinization is increasing globally in previously freshwater wetland environments (Herbert et al., 2015), causing large shifts in biogeochemical cycling as a result of salinity exposure (Herbert et al., 2015;Luo et al., 2019;Poffenbarger et al., 2011). The mechanisms behind salinization of tidal river systems and the adjacent terrestrial environments may be unique from the comparatively well-studied freshwater wetland salinization (Ensign & Noe, 2018). Salinization impacts not only biogeochemical cycling  but also coastal freshwater ecosystem function, including tree health . Understanding the mechanisms behind seawater inundation and floodplain salinization in coastal ecosystems is therefore necessary to elucidate how such chronic salinity exposure arises, persists, and is incorporated into predictive modeling efforts needed for projecting future change.
As land use can be a key driver of hydrobiogeochemical function in coastal ecosystems (Kroeger et al., 2017), land use alterations that impact hydrologic connectivity between tidal rivers and floodplains may also impact its biogeochemical functioning, potentially in similar ways as SLR. Anthropogenic changes to hydrologic flow such as dam building and removal represent a disturbance that has both short-and long-term impacts on hydrology and ecosystem function (Diefenderfer et al., 2008;Kroeger et al., 2017). For example, the increased frequency of inundation increases both sedimentation and channel areas in recently restored tidal wetlands (Diefenderfer et al., 2008) until the ecosystem reaches a new dynamic steady state.
While hydrologic connectivity across some common tidally influenced ecosystems, such as salt marshes, has been modeled with a variety of 1-, 2-, and 3-D models (e.g., see Table 1 in Xiao et al., 2017), the representation of terrestrial-aquatic fluxes of water and biogeochemical constituents is currently still missing from large-scale earth system modeling efforts (Clark et al., 2015;Fan et al., 2019;Pacific Northwest National Laboratory, 2019;U.S. DOE, 2017). Systematically mechanistic models that couple hydrology, biology, and chemistry are needed to understand coastal ecosystem function in terms of fluxes of water, carbon, and nutrients. In particular, incorporation of two-way interactions between terrestrial-aquatic systems, such as river-floodplain inundation and groundwater-surface water interactions are active areas of model improvement (Clark et al., 2015;Miguez-Macho & Fan, 2012a, 2012bYabusaki et al., 2017;Yamazaki et al., 2011). Here we observe and model the variably saturated flow and salinity in a floodplain of a first-order tidal stream, Beaver Creek, on the coast of Washington State, USA, to ascertain the hydrologic mechanisms behind the observed salinization of floodplain groundwaters following seawater flow restoration in the tidal stream in 2014. The watershed was restored to full tidal effects after nearly a century of freshwater conditions (Washington Department of Fish and Wildlife, 2019), which has had a marked influence on tree health  and soil organic carbon composition (Sengupta et al., 2019). We postulate that topography is a key determinant of floodplain inundation and groundwater salinization and that the two-way exchange of water in the tidal river system is dominated by surface connectivity and flood events.

Study Area
This study examines a tidal stream-floodplain system in western Washington State, USA, Beaver Creek. The stream is a 2.5 km first-order tributary draining into the larger Johns River (Figure 1), which drains into a nearby estuary, Grays Harbor, 1.5 km below the confluence with Beaver Creek. The Grays Harbor estuary has a~4.7 m tidal range (at NOAA Westport Station 9441102, inside the entrance to Grays Harbor, 10 km away from Beaver Creek). Salinity in the estuary varies seasonally from 1 to 33 practical salinity units (PSU) near the Johns River mouth (WA State Dept of Ecology, 2019), with highest salinities in late summer and lowest salinities in late winter. Grays Harbor is fed by several other, larger river systems, the most prominent being the Chehalis River (182 m 3 s −1 annual average discharge), 15 km northeast of the Johns River ( Figure 1).
A 1 m digital elevation model (DEM) of the region (LiDAR Bare Earth DEM, 2010) was used without correction to delineate the 382 ha Beaver Creek catchment ( Figure 1). Most of the catchment is classified as upland forest (Sengupta et al., 2019;Wang et al., 2019). The catchment consists of several lower-elevation floodplains with a consistent slope break around 3.7 m elevation (all elevations presented herein are NAVD88 unless specified otherwise) (Figure 1), accounting for~13% of the delineated watershed area (50 ha total). Beaver Creek is an ungauged watershed, therefore we estimated the annual freshwater discharge from (1) the annual rainfall (2,935 mm yr −1 ; US Climate Data, n.d.) at Hoquiam, (2) the delineated Beaver Creek catchment area (382 ha), and (3) the watershed efficiency (80%) from a gauged analog watershed, the Willapa River watershed.
In summer 2014, the Washington State Department of Fish and Wildlife removed a causeway to restore natural tidal movement of seawater to Beaver Creek for salmon and steelhead passage and rearing habitat (Washington Department of Fish and Wildlife, 2019). For nearly a century before causeway removal, Beaver Creek and its associated floodplains were a dominantly freshwater system. After causeway removal, saline water from the Grays Harbor estuary has been able to propagate through Beaver Creek during the upper limb of the rising tide. At low tide, freshwater baseflow continues to freely flow downstream. Floodplains near the mouth of Beaver Creek are now exposed to tidal water levels and salinities from the Grays Harbor estuary, whereas the headwaters and uppermost reaches are still freshwater. The focus of this study is on a floodplain of Beaver Creek that is between the freshwater and saline endmembers ( Figure 1). The floodplain is on the interior of a Beaver Creek meander bounded to the north by an upslope to higher elevations. The floodplain has a long axis of~250 m from southwest to northeast, and an average short axis of 75 m from southeast to northwest. The floodplain is relatively flat, with most of the central floodplain in the 2.5 to 2.8 m elevation interval. Higher elevations (>~3 m) can be found along the stream bank and the hillslope (Figure 2).

Water Resources Research
YABUSAKI ET AL.

Field Data
Water level and salinity were continuously monitored at one stream station and three floodplain wells along a transect from near the bank of Beaver Creek (18 m from the stream center), across the central floodplain (38 m from the stream center), to the hillslope (67 m from the stream center). Beaver Creek water level and water quality parameters were monitored at 30-min frequency from April 2018 to May 2019 using a multiparameter sonde (YSI EXO2, Xylem). Stream water levels were converted to common datum elevation accounting for dynamic fluid density and maximum water level across the transect during peak flooding conditions for each deployment segment. The near-bank, central, and near-hillslope floodplain wells were 1 m deep, 2-inch diameter PVC with 2.5 mm slots at 1 cm intervals across the entire belowground length (ESP Supply). Wells were screened using 250 μm well socks (ESP Supply) and sealed at the ground surface with bentonite clay to reduce artifacts of infiltration during flood events. Pressure and conductivity were monitored at 30-min intervals in each of the wells. Pressure head in the wells was determined from deployed pressure transducers~75 cm below ground level (HOBO U20L-04 Water Level, Onset Computer Corporation), correcting for barometric pressure at NOAA Westport Station (ID 9441102) at 30-min collection intervals from April 2018 to May 2019. Electrical conductivity readings from deployed devices at~60 cm below ground level (HOBO U24-002 Conductivity, Onset Computer Corporation) were converted to salinity following standard UNESCO/IOC/SCOR/IAPSO procedures as defined in the International Thermodynamic Equation Of Seawater-2010(IOC et al., 2010, using the Gibbs SeaWater Oceanographic Toolbox written for R ("gsw" package; Kelley et al., 2017). Importantly, this toolbox utilizes an updated Practical Salinity Scale-1978 (PSS-78) algorithm incorporating corrections for low salinities (Hill et al., 1986) and absolute salinity to describe the salinity of seawater (McDougall & Barker, 2011). Raw salinity was compensated for sensor drift and offset using device readings from seawater (Grays Harbor, WA, USA) and standard conductivity solution (50,000 μS cm −1 ; YSI) at beginning and end of each deployment. Salinity measurements presented in this work are represented by PSS-78 and presented in practical salinity units (PSU). Pressure head in wells was converted to hydraulic head accounting for fluid density and common datum elevation determined from flood events in July 2018, monitored by surface pressure transducers installed at ground level near each of the wells (HOBO U20L-04 Water Level, Onset Computer Corporation).
Soil samples at 10 and 30 cm (19 cm for the floodplain soils) were collected from the face of soil pits and soil classification (USDA Natural Resources Conservation Service, n.d.) and texture was determined. The floodplain soils are hydric Ocosta silty clay loam (USDA Natural Resources Conservation Service, n.d.). Soil texture was dominantly silty clay (n = 5), but ranged from sandy clay loam to clay (Sengupta et al., 2019).

Model Specifications
The objective of the numerical modeling was to mechanistically and systematically link observed behaviors in the targeted floodplain to the interplay between variably saturated flow and transport processes, floodplain material properties and topography, and dynamic site conditions. The PFLOTRAN simulator (Hammond et al., 2014) was used to model 3-D flow and salinity transport in the floodplain. A 3-D floodplain modeling domain was necessary to account for (1) low points along the stream bank that allowed inundation A single homogeneous material was specified for the entire model domain. Material properties (Table 1) included permeability, porosity, and van Genuchten water retention function parameters, which were guided by observed soil texture (Carsel & Parrish, 1988). The Mualem relative permeability function was used with the van Genuchten water retention parameters.
Flow in the floodplain was driven by Dirichlet pressure conditions along the stream and upland boundaries, and seepage face conditions on the ground surface. Stream boundary conditions were based on observed Beaver Creek stage whereas upland boundary conditions were based on observed water table elevations. Seepage face boundary conditions were used to represent stream interactions with the floodplain during inundation. When hydrostatic pressure head in the stream exceeded any ground surface elevation (i.e., top of ground surface grid cell), a ponding Dirichlet pressure boundary condition was applied; otherwise water was allowed to seep out of ground surface grid cells whenever internal water pressure exceeded atmospheric pressure. The two short segments on either side of the north boundary were treated as no-flow boundaries. Observed salinity in the stream was specified as a concentration accompanying water from the stream that entered the floodplain model domain via lateral flow across the aquifer-stream interface or infiltration through the ground surface during inundation events. The stream salinity was the only source of salinity to the modeled floodplain.

Model Spin-Up and Reconstruction of Floodplain Salinization
The model used the 1 year of observations from the stream and hillslope groundwater well as input. This hydrologic boundary condition data set was repeated for two model years to spin-up the subsurface flow model to a dynamic steady state of a repeating annual cycle of water levels, porewater saturation, and total water mass in the model domain. The resulting hydrologic initial condition was used to start the simulation of the multi-year salinization of the floodplain from an initially freshwater condition (i.e., prior to the restoration of seawater access to Beaver Creek in 2014). The salinization reconstruction comprises 4 years of continuous coupled flow and salinity transport between the creek and floodplain driven by repeated applications of the 1-year data set of boundary conditions to reach the period of observed data, which began on 23 April 2018. Longer-term floodplain salinization was simulated in a similar manner until 2038. Simulated floodplain inflows and outflows were continuously tracked to examine the water mass balance in the floodplain.

Observed Hydrology
Beaver Creek water levels were strongly tidal ( Figure 3). These water levels were highly correlated with measurements in Grays Harbor and display a 100-min time lag from the Westport gaging station, 10 km away from the stream station. Stream water levels represent the top half of the Grays Harbor tidal variation; that is, the stream bottom has a higher elevation than the lower 2 m of the tidal range in the estuary. Average freshwater discharge of Beaver Creek was estimated to be 0.28 m 3 s −1 . The average stream water elevation was 1.4 m (Figure 3a). Inundation events typically occurred with monthly regularity, according to the 29.53day lunar cycle. The magnitude and duration of floodplain inundation was seasonal, with the lowest monthly peak water levels (2.8 m) and shortest monthly inundation durations (9.5 hr) during summer ( Figure 3). The highest monthly peak water levels (3.3 m) and longest monthly inundation durations (44 hr) occur during the winter. Monthly mean water levels are higher in the winter and spring than in summer and fall, with the highest monthly mean water level occurring in February 2019 of the observed data set.
The floodplain water table elevation varied on a variety of time scales (e.g., seasonal, monthly, and tidal) yet generally maintained a water table gradient from hillslope to creek that was much steeper near the creek. The floodplain water table dynamics responded almost instantaneously to stream water levels greater than the local groundwater table, resulting in the water table near the bank to be most responsive to stream water levels ( Figure 3). The central floodplain well and the near-hillslope well only responded to progressively higher-inundation water levels (Figure 3). At the hillslope and bank locations, the monthly average water table was greater than 0.4 m below the ground surface. Inundation occurred only during the peak of the tidal cycle and lasted from 0 to 3 hr depending on the peak magnitude. The monthly tidal cycle typically included a few consecutive days with tidal inundation peaks ( Figure 3). Ponding depth and duration were highest where floodplain ground surface elevations were lowest, principally in the central floodplain. Monthly average floodplain water table elevations were seasonally variable; water table elevations were lowest in July and highest in December. This was in contrast to the stream, which was lowest in November and highest in February ( Figure 3).

Observed Salinity
Variability in Beaver Creek salinity was strongly correlated with tidal water level variability (Figures 3  and 4). At low tide during baseflow, stream water was essentially fresh and supplied from the upstream portion of the watershed. At high tide, salinities were maximal and reflected the variable salinity conditions of the greater estuary (WA State Dept of Ecology, 2019). Stream salinity exhibited a strong seasonality; with monthly average salinity peaking in summer for both the stream and floodplain groundwaters ( Figure 5). Stream salinities were highest in the summer when monthly peak water levels and inundation durations were lowest; stream salinities were lowest in the winter when peak water levels and inundation durations were highest (e.g., in August monthly peak creek water levels were 2.78 m and inundation duration was 9.5 hr vs. December monthly water levels of 3.31 m and inundation duration of 40.5 hr; Figure 5). Salinity peaks in the system were diminished by episodic rainfall events, which increased freshwater flow in Beaver Creek and watershed discharge into Grays Harbor as indicated by increased low tide water levels in the creek (Figure 3).
During floodplain inundation, stream water level and salinity was the highest of the tidal cycle, resulting in near-bank and central floodplain groundwater salinities that exceeded the average salinity in Beaver Creek (Figures 4 and 5). Floodplain groundwater salinity was variable both seasonally and spatially; lower salinities were observed with greater distance from the creek and during late fall through spring (Figures 4  and 5). Monthly average floodplain salinity was lowest near the hillslope, ranging from 1.6 PSU in May to 6.4 PSU in October 2018 ( Figure 5). The central floodplain monthly average salinity was highest for 8 months of the year, ranging from 6.1 PSU in May to 14.4 PSU in September ( Figure 5). The near-bank floodplain salinity was the most responsive to stream salinity, exhibiting a peak monthly average salinity of 15.0 PSU in August and minimum monthly average salinity of 6.0 PSU in December ( Figure 5). In November and December 2018, central floodplain groundwater salinities were higher than the average inundation stream salinity.

Floodplain Model
The 3-D model allowed for a detailed understanding of the importance of topography as a driver of floodplain salinity exposure. Model results suggest that during inundation, creek water fills the lower elevations in the central floodplain first, before expanding to higher elevations. Infiltration in the inundated areas is controlled by the antecedent moisture conditions, the proximity of the water table to the ground surface and ponding depth and duration, while groundwater salinity is controlled primarily by the infiltration and stream salinity during the inundation events. From the model, we observe a distinct topographical threshold of 2.6 m, where inundation of the central floodplain begins in advance of bank overtopping ( Figure 6). Furthermore, the spatially resolved model shows that although most of the floodplain is inundated after stream water levels reach 2.8 m, there continue to be substantial bank sections that are not submerged.

Large-Scale Total System Water Mass
Over the course of the year, the dynamics of the total water mass in the model domain is driven almost exclusively by infiltration during quasi-monthly inundation events during spring tides. The total water mass peaks and falls with the highest creek water level events, never recovering to a stable state or tendency toward a stable state of total water mass between events (Figure 7). Over seasonal time scales, total water mass trends lower from late spring through fall and higher from winter through spring, reflecting the higher extreme tide magnitude and frequency in the winter. Furthermore, on annual scales, total water mass variation shows no discernible interannual variation (Figures 7 and 9a), indicating that the system is in a repeatable dynamic steady state with respect to water flow. Over the annual cycle, the budget of water mass inputs Water Resources Research and outputs shows that over 97% of the water entering the floodplain comes from infiltration (average 3,660 Mg water over 3 years of simulation); less than 3% is derived from upland groundwater flow (average 110 Mg water over 3 years of simulation; Figure 7b).

Modeled Floodplain Salinity
The modeled floodplain is assumed to have been a freshwater ecosystem prior to restoration of seawater access to Beaver Creek in 2014. This is based on an observed rapid decline in tree growth rates immediately following culvert removal . The observed data set from April 2018 to May 2019 is assumed to be Year 4 of floodplain salinization. Subsurface floodplain salinity is laterally controlled by topography, as illustrated by the spatial salinization pattern of the floodplain ( Figure S2). The subsurface floodplain salinity is also vertically controlled by ponding and antecedent moisture conditions and permeability of floodplain  Total salt mass dynamics are also driven by inundation events, similar to water mass. The seasonal trend is opposite from total water mass; total salt mass increases from spring through early fall, declines from late fall

10.1029/2019WR026850
Water Resources Research to early winter, and stabilizes through early spring at a higher level than the previous year. Furthermore, there is considerable salt accumulation from year to year, and after the 4-year reconstruction of floodplain salinization, salt mass is still accumulating (Figure 9b). The total salt mass variation reaches a near-dynamic steady state condition after 20 years of simulation, when the salt efflux to the stream begins to balance the continual salt additions (Figure 9c). This is essentially the transport time for salt entering the floodplain near the hillslope to finally begin reaching the creek.

Data-Model Comparison
Model goodness-of-fit was evaluated using a variety of common evaluation statistics, group statistics and graphical evaluation techniques (Legates & McCabe, 1999;Moriasi et al., 2012Moriasi et al., , 2015Nash & Sutcliffe, 1970;Willmott et al., 2015), based on observations at the three floodplain well locations. The floodplain model is able to capture many characteristics of observed floodplain well water levels and salinities across the floodplain (Figures 10, 11, and S3). An evaluation of model fit using the Willmott's refined index of agreement (d r ) (Willmott et al., 2012(Willmott et al., , 2015 indicates that over the year, the model-estimated variability is an acceptable estimate of the observed variability for both water level and salinity across the floodplain (Tables 2 and 3; Figures S4 and S5). The yearly d r ranged from 0.39-0.67 for water level and 0.24-0.63 for salinity, where −1 represents the highest mean absolute error relative to the mean absolute deviation, and 1 represents the lowest (Tables 2 and 3). Additional yearly model evaluation metrics examined herein are available in the supporting information (Tables S1 and S2). Notably, evaluation of model performance was highly

10.1029/2019WR026850
Water Resources Research dependent on the statistical metric applied, highlighting the sensitivity of commonly used model evaluation statistical metrics to highly variable data sets (Schaefli & Gupta, 2007;Willmott et al., 2012), and the importance of direct data-model comparisons for such data sets (Bennett et al., 2013;Moriasi et al., 2012).
The model is able to replicate the water table profile from near-hillslope through the central floodplain to the near-bank; the timing of water level dynamics over semidiurnal, lunar monthly, and seasonal time scales; the magnitudes of the highest water level peaks; and the water table behavior between inundation events (Figures 10 and S3). The yearly mean, median, and standard deviation of modeled results closely follow the observations across the floodplain (Table 2), highlighting overall model performance on yearly time scales. Hourly, semidiurnal, diurnal, and monthly mean absolute errors are similar, indicating model validity over these time scales (Table S3). The model is particularly effective at reproducing the observed winter water level dynamics when inundation events are the largest and most frequent throughout the floodplain, highlighted by lowest mean error during those tidal cycles ( Figure S4). The most noticeable discrepancy is an underprediction of the near-bank water levels between summer inundation events (Figure 10), where the model performance as assessed through d r is the worst during the new moon cycles during this period ( Figure S4). At this location, the model is less influenced by the secondary water level peaks between the larger monthly peaks than the observed data.
The response of the modeled floodplain water levels to the boundary conditions is largely dictated by the permeability specification of 4.86 × 10 −13 m 2 (hydraulic conductivity of 3.65 × 10 −6 m s −1 ), which controls the unconfined water table gradient and infiltration flux during inundation. This permeability specification systematically honors the observed water levels at multiple floodplain locations over a range of important time scales (e.g., semidiurnal tides, monthly lunar cycle of inundation, seasonal variation) and thus allows for model assessment of floodplain behaviors that are not measured (e.g., flow, velocity, infiltration, saturation, total water mass).
The model captured seasonality of salinity observations across the floodplain, in addition to magnitude in the near-bank and near-hillslope floodplain wells (Figures 11, S3, and Table S3). However, for the central floodplain well the model generally underpredicts the salinity magnitude by 4 to 8 PSU. This results in the highest yearly mean absolute error of 3.79 PSU (Table 3). Similar to modeled floodplain water levels, the model salinities tend to respond less dynamically to the inundation events than the well observations, resulting in a range of model performance across tidal cycles and locations ( Figure S5). Despite this, annual mean, median, and standard deviation of the modeled salinity is well matched to observed salinity (Table 3), with similar mean absolute error over hourly, semidiurnal, diurnal, and monthly (Table S3) scales, confirming model validity over a variety of time scales.

Parameter Sensitivity and Uncertainty
Simulations were performed to test the sensitivity of the simulated floodplain water levels and salinity to the model specification. The prior uncertainty associated with the four parameters are fully explored using an efficient sampling-based uncertainty quantification (UQ) framework. Uniform distributions are assumed to

Water Resources Research
consider the full ranges of prior uncertainties. The four-dimensional parameter space was sampled using the Latin Hypercube Sampling (LHS) method, resulting in the generation of 64 combinations of samples (i.e., realizations). The advantage of LHS over Monte Carlo sampling is that the samples can fully cover the parameter space with even spread, minimal gapping, and minimal lumping (i.e., duplications) to achieve both high efficiency and representativeness. This is illustrated in scatter plots of the sample sets shown in Figure S6.
Generalized linear model (GLM) fitting and analysis of variance (ANOVA) for the sample combinations showed that~80% of the total variability in hourly water levels and~60% of the total variability in hourly salinity can be explained by the GLM models with all four parameters included (Table S4). Among the explained water level and salinity variability, log permeability alone is contributing to~90% and~80%, respectively. The secondary but nonnegligible parameters are porosity and stream water level boundary condition offset, while the hillslope groundwater level boundary condition offset seems to be a trivial factor for the current list of response variables. In this assessment, the parameter importance ranking can guide future data acquisition to maximize the information gain toward model improvement, and the ensemble misfit analyses provide ranges of parameter values corresponding to low misfit, which can be useful for calibrating the mechanistic models toward higher accuracy.
Despite being a dominant parameter, there was no single value of permeability able to reproduce in detail the behaviors at all floodplain observation well locations: near bank, central, and near hillslope. In general, the high end of the hydraulic conductivity range tends to underpredict water table elevations between peaks at all three locations; whereas the predicted water tables using the low end of the hydraulic conductivity range tend to match the near-bank floodplain water levels, while overestimating the central and near-hillslope water levels. Similarly, there was no systematic improvement of model fidelity with any combination of water elevation offsets. Typically, a given offset combination would improve fidelity at one floodplain well location at the expense of fidelity at one or both of the other floodplain well locations.
Accurate topography is very important to floodplain salinity infiltration. The ground surface elevation determines inundation occurrence, ponding depth and duration, and depth to water table. This can be important when comparing model results with field observations. In this sensitivity analysis, results from grid cells near the floodplain sampling well locations but with different ground surface elevations were compared to the grid cell containing the well. Near the hillslope well, a 0.3 m lower ground surface elevation had as much as a 10 PSU higher salinity. Near the central floodplain well, a 0.1 m higher ground surface elevation had as much as a 4 PSU higher salinity. While these results highlight the role of topography on salinization dynamics, no nearby grid cells had results that were substantially closer to observations than the cells where the actual wells were located. Future work with increased model elevation resolution may help further elucidate the impact of microtopography on salinization dynamics.
The impact of salinity density effects was tested on a 2-D vertical cross section from the upland boundary to the stream along the transect of the floodplain wells. Floodplain water levels and salinity, driven by the dynamic stream and hillslope boundary conditions, were simulated with and without salinity density effects invoked. The resulting spatiotemporal salinity distributions were indistinguishable. Consequently, salinity density effects were omitted from the simulations generating the presented results.

Drivers of Terrestrial-Aquatic Connectivity in Tidal Streams
In the Beaver Creek watershed, floodplain inundation is topographically controlled, similar to many other observed coastal ecosystems (Xin et al., 2013). The modeled total system water mass responds primarily to inundation events, similar to porewater turnover in tidal marshes, which are dominated by estuarine surface water exchange (Harvey & Odum, 1990). This similarity is consistent with the evolution of the Beaver Creek floodplain from a freshwater swamp to a marsh ecosystem as the floodplain becomes saline . During the largest infiltration events, the mass of infiltrating water is less than 1% of the overlying inundation water mass. Thus, most of the inundation volume drains back into the creek as overland flow, rather than completely infiltrating into the floodplain soils and laterally flushing across the aquifer-stream interface. This is in contrast to larger systems such as the Yangtze River and Chesapeake 10.1029/2019WR026850

Water Resources Research
Bay, where surface seawater infiltration is primarily discharged via groundwater pathways along the tidal river bank, creating a dynamic infiltration-flushing circulation of seawater (Harvey et al., 1987;Xia & Li, 2012;Xiao et al., 2017), and the modeled results of increased groundwater flow from macrotidal systems in nonriver dominated estuaries with SLR (Wilson & Gardner, 2006).
The flow of water through the terrestrial-aquatic interface is largely controlled by soil permeability/hydraulic conductivity, which is a primary determinant of surface water-groundwater exchange rates in coastal systems (Hemond & Fifield, 1982;Moffett et al., 2012). This importance is reinforced by the substantial effect of permeability on simulated water levels and salinities identified in the UQ analysis. In our modeled system, similar to many coastal floodplains, floodplain soil permeability is low, limiting groundwater-surface water interactions and lateral transport (Harvey et al., 1987). The permeability is sufficiently low to limit lateral penetration of stream water across the aquifer interface, despite the~2.5-m range in stream water levels.
Low permeability also limits infiltration during inundation events to a small fraction of the overlying inundation volume; however, based on the total mass dynamics, infiltration is the principal interaction between the stream and floodplain groundwater, as well as the mechanism for salinity delivery to the floodplain at large. As a result, the primary dynamics in the floodplain are vertical as exemplified by subsurface salinity distribution. Heterogeneity of actual soil hydrologic conductivity across a system (Gupta et al., 2006) and/or fast pathways, such as surface and subterranean erosion channels (Dunne, 1990) or macropores (Johnston et al., 2005;Xin et al., 2016), are plausible explanations for observed dynamic response of water levels across the floodplain and may explain the observed lag in the model water table response, warranting further examination. However, while potentially important, these processes are likely not dominant in this system as we can largely reconstruct the observed data with our current model efforts. Further considerations of such microscale heterogeneity could further validate and constrain the model with additional data assimilation of variables such as spatially and depth-resolved soil moisture/conductivity sensors, spatial and depth-resolved soil properties and finer-resolution topographical representation.
The apparent drivers of inundation dynamics highlight the temporal and spatial variability of floodplain inundation, creating hot spots and hot moments of seawater exposure in the restored floodplain. Spatially, inundation dynamics may be related to channel-like morphology, such that smaller channel dynamics may be a major mechanism for tidal hydrologic flow and groundwater drainage in this system, similar to tidal marshes (Moffett et al., 2012;Xin et al., 2011). In fact, bank overtopping is not necessary for inundation across the floodplain. Temporally, the mechanisms behind the observed seasonality of water infiltration may also be linked to a balance of hydrologic processes, such as the interplay between shifts in rainfall, hillslope inflow and flood frequency driving water table dynamics (Burt & Pinay, 2005) or the interplay between evapotranspiration and tidal forcing (Hemond & Fifield, 1982). For example, inundation of the floodplain varies not only on spring/neap tidal time scales but also seasonally. In the summer, the monthly average floodplain water table is at a seasonal low, due to the smallest inundation peaks and shortest inundation periods of the year. The lower water table provides capacitance for the infiltration of stream water during inundation events, as infiltration can occur once stream surface water level exceeds groundwater table elevation in the surrounding floodplain (Harvey et al., 1987). In the fall, the monthly average creek water level is similar to the summer water levels, but the duration of floodplain inundation increases. This leads to increasing water table elevations that plateau in the winter. From late fall through spring, the higher floodplain water table and ground saturation limits further infiltration from subsequent events. Interestingly, the general hydraulic gradient from hillslope-floodplain-stream is retained year-round, despite reversal of hydraulic gradients and disconnection across landscape units observed in nontidal floodplain systems (Burt & Pinay, 2005) and zero groundwater head gradient in coastal marsh ecosystems (Nuttle, 1988).

Floodplain Salinization
The infiltration dynamics are the dominant driver of floodplain salinization at Beaver Creek. The temporal nature in stream salinity coupled with inundation magnitude, frequency, and duration result in a dynamic salinization of the floodplain. Depth to water table in the floodplain can be small and the vadose zone is often saturated or nearly saturated during inundation events. This limits infiltration and explains the observed persistence of central floodplain salinity during the highest and longest low-salinity inundation events of the winter. Therefore, while stream salinity during inundation has a 10.1029/2019WR026850

Water Resources Research
prominent impact on floodplain groundwater salinization, its immediate influence on groundwater salinity is regulated by soil permeability, ponding depth, duration, and antecedent moisture conditions. The low permeability of floodplain soils in our system results in low lateral subsurface flow and long groundwater residence times, which has ramifications for the post-restoration accumulation of salt in the floodplain groundwater. Furthermore, these soil properties suggest that salinization mechanisms of this ecosystem are more similar to that of marsh systems than that of coastal aquifers/groundwaters. For example, compacted subsoils or confining layers severely limit groundwater salinization (Crooks & Pye, 2000;Schultz & Ruppel, 2002).
Based on our model, the time required for the annual loading of salt to the floodplain to equal the annual efflux of salt to the river is over 20 years from the restoration event (i.e., barrier removal in 2014). Saltwater intrusion into floodplain groundwaters thus compound press and pulse salinization dynamics. Press salinization in the system is resultant of the year to year accumulation of salt in the floodplain soils, revealed by the salt mass balance. Pulse salinization is based on seasonal and spatial dynamics of observed and modeled salinity; along the transect from the hillslope to the creek, the vertical salinity distribution is localized to lower floodplain elevations. The amount of time salt mass in the floodplain is greater than year-end net accumulation increases throughout the simulated period from about 30% of the year in 2018 to over half of the year in 2038, indicating that pulse events remain a persistent and increasingly important component of salinization in this system. Both press and pulse salinization affect ecosystem processes and biogeochemical cycling, although perhaps on differential scales (Chambers et al., 2013;Herbert et al., 2018;Weston et al., 2011). For example, fluctuating salinity incursions in coastal freshwater wetlands affect redox regimes and sediment biogeochemistry, while salinity exposure impacts microbial utilization of organic matter (Schoepfer et al., 2014;Sengupta et al., 2019;Weston et al., 2006). The observed and modeled mechanisms of salinization in this ecosystem may have additive or differential effects on floodplain biogeochemistry (e.g., Neubauer, 2013) and highlight the importance of monitoring and elucidating compounding drivers of ecosystem process and function.

Restoration of Tidal Flow Foreshadows the Impact of Natural Disturbances
On yearly time scales, dynamic flow and salinity behavior in Beaver Creek are driven by vertical flow and tidal impact, while on decadal time scales, lateral transport across the floodplain is necessary to reach dynamic steady state for total salt mass. The post-restoration response of the system represents a press response to a change in hydrologic regime with long-term impact on the floodplain salinity. The influence that restoration events have on salinization of floodplain groundwaters are not the same across systems and may be driven by groundwater hydraulic properties (Johnston et al., 2005); groundwater exchange and salinization are varied in response to restored tidal access and appear to be impacted by the system geomorphology (e.g., macropore presence, soil hydraulic conductivity, soil structure; Johnston et al., 2005;Van Putte et al., 2019). Thus, restored sites such as Beaver Creek may provide useful analogs to mechanistically understand the response of coastal freshwater ecosystems to increasing inundation and salinity from press and pulse disturbances, such as SLR, storm surge events, land use change, and management practices.
Understanding the time scales of coastal ecosystem salinization is critical for determining the potential impact of these numerous drivers of inundation regime shifts on both coastal ecology and biogeochemical function. Coastal forests around the Pacific Northwest, for example, have been consistently declining in health with elevated tree mortality and declining growth rates over the last several decades due to press salinization, while events such as the Beaver Creek culvert removal can have the same effects in a year or less time . Similar observations of long-term coastal forest mortality due to SLR have been made on the U.S. East Coast (Kirwan & Gedan, 2019), and likewise, extreme events such as hurricanes can have similar effects as to what we observed after culvert removal (Cahoon et al., 2003). These long-term and rapid changes in ecology can have a significant impact on biogeochemical cycles. For example, tree mortality after Hurricane Katrina on the U.S. Gulf Coast damaged a reservoir of carbon equivalent to 50% to 140% of the net annual U.S. forest carbon sink (Chambers et al., 2007). Changes in hydrological dynamics underscore all the diverse disturbances described above. Thus, our ability to mechanistically predict the nature of anthropogenic and natural changes in hydrological regimes is paramount in understanding how Earth systems will respond from local to global scales.

Conclusion
Systematic and mechanistic numerical modeling generated new understanding from a 1-year data set of the first-order tidal stream-floodplain system. While it may have been possible to use the observations to qualitatively infer the link between floodplain groundwater and monthly inundation events, the integration of modeling and measurement enabled the quantitative assessment of many important features that were not directly measured. Foremost was the conclusion that the Beaver Creek floodplain is in an early stage of salinization and that salt will continue to accumulate for~20 years. This was due to the mechanistic interplay of topography, low-permeability materials, soil physics, and the seasonality of the stream water levels and salinity. While total water and salt mass clearly show the controlling impact of monthly inundation events, the topography and local hydrologic conditions regulate the spatial sensitivity to those events. This results in considerable spatial variation in salinization across the floodplain.
These and other emergent behaviors determined by the modeling underscore the value of Beaver Creek and other freshwater coastal sites in transition after restoration of seawater access. These sites can be useful natural analogs to study the impact of encroaching seawater exposure on coastal freshwater ecosystems in anticipation of SLR, increasing storm surge, climate change, and land use change. The decadal transition time scale provides a research opportunity to observe biogeochemical response to inundation and salinization in real time. Furthermore, the spatial variation in floodplain response can be used in a space-for-time analysis of early stage impacts of seawater exposure.

Data Availability Statement
All data are provided in the figures and tables of this manuscript; full data from the model outputs can be found at the figshare.com repository (10.6084/m9.figshare.c.4758797.v2). Model code is available in the cited references and at this site (https://www.pflotran.org/documentation).