Forecasting Sea Level Rise-driven Inundation in Diked and Tidally Restricted Coastal Lowlands

Diked and drained coastal lowlands rely on hydraulic and protective infrastructure that may not function as designed in areas with relative sea-level rise. The slow and incremental loss of the hydraulic conditions required for a well-drained system make it difficult to identify if and when the flow structures no longer discharge enough water, especially in tidal settings where two-way flows occur through the dike. We developed and applied a hydraulic mass-balance model to quantify how water levels in the diked and tidally restricted coastal wetlands and water bodies dynamically respond to sea-level rise, specifically applied to the Herring River Estuary in MA, USA, from 2020 to 2100. Sensitivity testing of the model parameters indicated that primary outcomes were not sensitive to many of the chosen input values, though the terrestrial water input rate to the estuary and the flow coefficient for the hydraulic infrastructure were important. The relative importance of parameters, however, is expected to be site specific. We introduced a drainability metric that quantifies the net water volume drained over every tidal cycle to monitor and forecast how rising water levels on either side of the dike affected the net draining or impounding conditions of the system. Ensembles of model results across parameter and sea-level scenario uncertainties indicated that substantial impoundment of the Herring River Estuary was expected within ~ 20 years with the existing flow structures, a sluice and two flap gates. Simulations with up to three additional gates did not dampen this trend toward impoundment, suggesting that rising impounded water levels are likely even with major construction upgrades. Increasingly impounded diked coastal waterbodies present a hydrologic challenge with socioecological implications due to projected flooding and ecosystem impacts. Solutions to this challenge may be to allow coastal wetland restoration pathways or require substantial and recurring infrastructure improvement projects.


Introduction
A substantial fraction of coastal wetlands today is impacted to varying degrees by human activities (Gedan et al. 2009;Burdick and Roman 2012;Pendleton et al. 2012;Kroeger et al. 2017;Crooks et al. 2018). In both urban and rural settings, wetlands that have not been destroyed by filling or dredging commonly persist within disturbed hydrologic regimes. Such disturbances result from diking to restrict tidal inundation, ditching to drain either tidal or meteoric water more efficiently, and other alterations. Coastal hydraulic structures, such as dikes, are designed for the hydrologic setting when built. Performance of such infrastructure may decline as climate change, subsidence, and sea-level rise alter the hydrologic regime and raise the tidal frame (Spencer et al. 2016;Schuerch et al. 2018). Due to sensitivity to hydrology, changes in hydrologic regimes are likely to drive changes in the ecological and biogeochemical functioning of diked ecosystems, creating positive feedbacks in both climatic and hydrologic drivers by altering greenhouse gas fluxes and relative soil elevation (Portnoy 1999;Wang et al. 2019;Eagle et al. 2022). Further, aging dikes designed for historic Communicated by David K. Ralston hydrology could exacerbate flooding hazards and damage coastal infrastructure. Predictions and models are lacking for the hydrology, ecology, and carbon cycle processes in diked landscapes, and thus, as sea levels rise, the future hydrologic states and fates of widespread coastal diked lands are critical, and yet highly uncertain.
Coastal dikes are intended to protect low-lying areas from inundation or short-term flooding from marine, often tidal, waters. However, dikes can also impound freshwater flows and tidally exchanged water caused by storms and/or wave overtopping, requiring adequately sized outlet flow control structures to drain stored water. Depending on the diked system, the flow control structures can be limited to draining only during low tides with one-way features, such as flap gates, effectively maintaining low salinity in impounded water or discharge saline water intruded during extreme events. If some tidal exchange is allowed across the dike, culverts or sluices (i.e., gate with an adjustable height opening at the bottom) can be sized to limit the volumes of exchange. One-way flow control structures also commonly lead to more stable and overall lower water levels on the inland side of the dike and a drier landscape. We refer to these two inland diked conditions as either impounded (i.e., insufficient outflow) or drained (i.e., enhanced outflow and restricted inflow) dike systems, depending on the degree to which diked water levels remain consistently higher or lower than a natural tidal setting.
We anticipate that most dike systems were designed to protect or alter a specific area, such as by draining upland areas and restricting tidal inundation to expand coastal agricultural land and reduce mosquito habitat (Kroeger et al. 2017;Crooks et al. 2018). Higher sea levels would reduce the amount of time during each low tide for diked systems to drain. Among our hypotheses is that, without changes to the dike infrastructure, such reductions in drainage will eventually convert drained systems into impounded systems. This transition may take decades to manifest and may be obscured by tidal and sea level dynamics. Watershed morphology, hydrology, weather, and infrastructure complicate the magnitude and timing of the hydrologic changes and combine to create unique conditions for how each diked system will respond.
Transitioning hydrologic regimes feedback into coastal ecosystem dynamics, blue carbon cycling, and associated greenhouse gas emissions. With additional freshwater impoundment, increased terrestrial inundation and more extensive freshwater features could replace brackish or saltwater wetlands, expanding the potential for methanogenic conditions (Warren et al. 2002;Kroeger et al. 2017;Sanders-DeMott et al. 2022). In addition, while drained former coastal wetlands emit sequestered soil carbon, it is not straightforward to predict the impact of a transition to an impounded condition on the rate of carbon storage, since impounded coastal wetlands have suppressed rates of soil carbon storage relative to rates in natural tidal wetlands (Turner 2004;Eagle et al. 2022).
The focus of this study is to investigate how coastal diked systems respond to sea-level rise and test how the drainage performance of existing hydraulic infrastructure changes in the future. We develop a simplified water balance approach to model flow through hydraulic control structures and water level in a diked system. We test this model with observations from the Herring River Estuary in MA, USA, for presentday and future sea-level rise scenarios. We hypothesize that drained diked systems will become impounded with sealevel rise if existing infrastructure remains functional but not upgraded. For our study site, we forecast when the system converts from drained to impounded, and we test a methodology to perform such an analysis on other diked systems. We further demonstrate the potential of several alternative water control structure scenarios to change the impoundment forecasts of the Herring River Estuary.

Site Description
Our study used the diked Herring River Estuary (HRE; Fig. 1) watershed to develop and test the model framework. The HRE watershed spans ~ 19 km 2 of sandy glacial till deposits. A dike was constructed across the HRE in 1909, resulting in reduced tidal inflow, lower average water level, and drainage of the HRE (Portnoy and Giblin 1997;Portnoy 1999). The current flow controls consist of two identical flap gates and one static sluice gate that was historically adjustable. The flap gates allow the HRE to discharge to Wellfleet Harbor at lower tidal conditions and are otherwise closed. Wellfleet Harbor tides at the dike are semidiurnal and mesotidal with a principal lunar semidiurnal (M2) amplitude of 1.23 m (Mullaney et al. 2020). The elevation of the top of the dike is 3.6 m above the North American Vertical Datum of 1988 (NAVD88), which was also used as the datum for all other elevations and water levels. The sluice allows flow in both directions resulting in a dampened tidal response in the HRE with an M2 amplitude of 0.33 m (Mullaney et al. 2020). The underlying unconfined aquifer discharges an annual average of approximately 0.03 m 3 /s of baseflow to the HRE based on numerical simulations (Masterson 2004).
Water levels have been recorded at a 5-min interval on each side of the flow control structure from September 2017 to present (site number 011058798)(U.S. Geological Survey 2022). Measurements of water level and discharge began in June 2015, but the discharge measurements were discontinued in September 2017 when a tidal radar water level measurement device was installed on the downstream side of the flow control structure. Without coincident discharge observations with water level measurements on each side of the dike, we only tested the model performance with the water levels and not with discharge. Note, however, that the water level measurements on the harborside of the dike are raised by the discharge from the HRE during low tide, such that the water levels at this location remain above low tide levels farther offshore from the dike within Wellfleet Harbor. This water level difference with the harbor required an approximation for forecasting water levels using tide gauge measurements and sea-level rise rates for Cape Cod Bay described in the "Sea-Level Rise Forecasts" section.

Methods
In the following sections, we develop the formulation of the hydraulic flow model within the context of tidally restricted and diked coastal systems ("Hydraulic Flow Model" section), we establish how the model results for hindcast simulations are tested for accuracy while quantifying the uncertainty of the outputs using Sobol sensitivity testing ("Sobol Sensitivity Testing" section), and we introduce a framework for applying and analyzing the hydraulic flow model with future sea-level rise and alternative infrastructure for the HRE ("Sea-Level Rise Forecasts" section).

Hydraulic Flow Model
We developed a hydraulic mass-balance model for the flow through the two flap gates and sluice. The objective of the model was to solve for the water level in the HRE with the water level in Wellfleet Harbor as the main boundary condition. For the development of the model, we selected a present-day 2-year period with continuous water level observations on either side of the dike control structure, from 27 February 2020 to 27 February 2022. First, water level observations were downsampled from a 5-min interval to a 30-min interval, balancing a high enough resolution to simulate tidal levels and flow with computational efficiency. Second, we modeled the volumetric flow rate through the control structure using a series of hydraulic equations based on the relative water level difference between the HRE and the outflow to Wellfleet Harbor. Finally, the water level in the HRE was updated for each 30-min time step using a stage-storage relationship for the HRE upstream to the embankment of High Toss Road, up to approximately the 5-m elevation contour, and the dike (Fig. 1). The stage-storage relationship was developed with a low-tide 1-m spatial resolution elevation model with additional surveyed bathymetry extracted from a hydrodynamic model (WHG Inc. 2012) and higher elevations bilinearly interpolated from a 10-m resolution digital elevation model (U.S. Geological Survey 2018). The stage-storage relationship was developed across the elevation range − 2.5 to 3.6 m with the low value set below the invert elevation and the latter equaling the height of the dike. The net flow into or out of the HRE was calculated for each time step by summing the flow through each dike control structure with the amount of streamflow and groundwater discharge entering the HRE from its contributing watershed, q in . We kept q in constant in time, but it was varied between model runs, which was estimated in an earlier analysis to be a constant 0.34 m 3 /s (WHG Inc. 2012). Once flow through each culvert bay was calculated, the volume of water stored in the HRE was calculated by multiplying the net flow rate by the time step. Finally, the stage-storage relationship provided the new HRE water level for the calculated volume.
We modeled the sluice gate using a series of hydraulic equations for five flow regimes that could occur in either flow direction (Fig. 2). The relative water levels on either side of the sluice gate controlled the flow direction and the flow regime with the relationships provided in Fig. 2 based on the implementation of sluice and weir flow in the HEC-RAS River Analysis System (Brunner 2016). All the equations for the flow followed the form (Brunner 2016): where Q is the volumetric flow rate of water flowing across the sluice, C is a unitless discharge coefficient, A is the cross-sectional area of the sluice or gate opening, H is the relative height of water driving flow, n is a unitless flow factor based on the hydraulic setting equal to one unless otherwise noted in Fig. 2 (Brunner 2016), and g is Earth's gravitational acceleration. For low water levels in the sluice, A is dependent on H. When submerged, A was calculated by the multiplication of the sluice gate width, w = 1.829 m, and its opening height, B = 0.485 m. For lower harborside water levels, A is equal to H multiplied by w. The height of the invert above the harborside bed, P, was not included in the calculations (Fig. 2a) where H loss decreases linearly from the maximum head loss, H loss_max , as the water depth on the marine, Z M , and land side, Z L , approaches a minimum water depth difference, D HL , where head losses created by the flow become negligible. Head losses in the culvert leading up to the gates were not considered. The two rectangular flap gates were modeled using a moment balance to calculate the flap opening angle, θ, created by water flowing out of the HRE acting against both the weight of the flap and hydrostatic pressure on the outer side of the submerged portion of the flap with: The moment balance included the weight of the gate materials, W, assumed to be uniformly distributed over the gate having a vertical length of h flap . The forces of the water acting on each side of the gate were set by the wetted width of the gate on the Wellfleet Harbor side, w out , and on the HRE side, w in . On the Wellfleet Harbor side, the height of the gate hinge above the harbor water level at a given time was d down . On the HRE side of the gate, the height of the gate hinge above the estuarine water level at a given time was d up . The flap gates were not vertical when closed, having a starting angle, θ 0 . We used the density of seawater, ρ, for water flowing out of the HRE, as the sluice gate allows seawater inflow during high tides that maintains high salinity across tides from salinity observations at the dike. The flap opening angle was solved numerically by finding the root of the zero-sum moment balance equation with the opensource Python SciPy package nonlinear function, optimize.fsolve (Virtanen et al. 2020). Once the flap opening angle was solved, the flow through the flap opening was modeled with weir equations (Eq. (1)) with an additional headloss term in H loss (Eq. (2)). No flow was allowed to occur through the flap gates when they were closed (i.e., θ = 0), although some leakage through the aging flaps was observed during a site visit in 2019. (2)

Sobol Sensitivity Testing
A global sensitivity test was performed for the HRE application of the flow control structures with the Sensitivity Analysis Library in Python (Herman and Usher 2017). Flap gate and sluice flow coefficients, head loss terms, initial estuarine water level, and q in were sampled for individual Sobol sensitivity tests using the extended form of the quasi-random samples of the parameter space including second order sensitivities (Saltelli 2002). The sampling routine was used 2 10 times with the 18 hydraulic and initial condition parameters to produce 38,912 parameter combinations. Each parameter was allowed to range between a value of 0.01 to 2 with parameter-specific units. The values and ranges of C are dependent on the type, construction dimensions, and hydraulic regime (i.e., conditions in Fig. 2), where C for sluice gates is often in the range of 0.5-0.7 and for weirs between 2.6 and 4.1 (Brunner 2016). For infrastructure with available specification sheets, the C value for certain flow conditions may be available for different construction geometries from experiments performed by the manufacturer. Previous work at the HRE found that C under the weir conditions remained below 2.0 (WHG Inc. 2012), but testing with higher C values may be needed in other systems. Similarly, q in is dependent on the site hydrology, climate, and size of the watershed, where an upper limit of 2 m 3 /s would only apply for relatively small systems like the HRE. Since q in integrates the combination of all rivers, streams, and groundwater to the diked or tidally restricted system being modeled, the range of q in required will depend on the total freshwater influx. For each parameter combination, the observed water levels within the HRE and outside of the dike in Wellfleet Harbor were run through the model with a 10-min time step over 30 days to calculate a non-parametric model efficiency, R NP (Pool et al. 2018): with: β is the ratio of the mean observed, μ obs , and mean simulated, μ sim , water levels in the HRE. The normalized variability, α NP , in the HRE simulated and observed water levels, h sim and h obs , used the difference between the ranked water levels, I and J. The Spearman rank correlation, r s , described the similarity in the dynamics of the water levels using the difference between the rank, R sim and R obs , and the mean rank, R sim and R obs , in the simulated and observed time series, respectively. The sensitivity of the efficiency criteria to the flow coefficients can be used to identify which parameters influenced the accuracy of the model. Multiple output parameter combinations that yield high efficiency values (e.g., R NP > 0.8) demonstrate the non-uniqueness of the hydraulic parameters in the model. These parameter combinations could be used for the remainder of the analysis for a quantification of model uncertainties. We adopted the R NP as the goodness-of-fit criteria, although other efficiency metrics could serve the same overall purpose to constrain hydraulic parameter uncertainties. In addition to R NP , we calculated the Nash-Sutcliffe efficiency, the relative efficiency, and the Kling-Gupta efficiency during model development with little effect on the total number of "good" model results (Nash and Sutcliffe 1970;Krause et al. 2005;Gupta et al. 2009). We focus on how these parameter combinations are tied to the implications of sea level changes on the operations of the existing and alternative flow control structures. Hence, we use multiple parameter combinations with high efficiency criteria from the sensitivity testing, with "good" defined as R NP > 0.9, in the following analyses to constrain uncertainty introduced by the model parameters.
We further tested these parameter combinations for the present-day 2-year period with water level observations, 27 February 2020 to 27 February 2022, at the longer 30-min time step used for the forecasting models. We calculated a new R NP for the longer model for all R NP > 0.9 models from the sensitivity testing. Models that performed well, with R NP > 0.8, for the 2-year simulation were included in the sea-level rise forecast models.

Sea-level Rise Forecasts
To quantify how the HRE water levels respond to sea-level rise, we constructed synthetic water level time series on the Wellfleet Harbor side, hereafter referred to as harborside, of the flow control structure. To accomplish this, we first calculated the amplitudes and phases of the tidal harmonic (7) constituents with a least squares regression from the dikeadjacent harborside water level observations during the 2-year present-day simulation period (27 February 2020 to 27 February 2022) with the Python package pytides (https:// github. com/ sam-cox/ pytid es). Because the harborside water levels are influenced by discharge through the dike, the tidal reconstruction created lower low tide water levels compared to observations (Fig. 3a). We developed a method to restrict low harborside water levels based on an apparent linear relationship between the minimum harborside and HRE water levels. This rough linear relationship physically represents the role of relative water level differences controlling discharge to the harbor and setting the water level of the outflow. Thus, if both water levels are known, the harborside water level could be maintained above a minimum elevation and used to correct harborside tidal reconstructions. For the present-day simulation period, we found little change in the model performance for the tidal reconstructions with and without the minimum harborside water level linear relationship. Given this similarity in present-day results and the need for HRE water levels to be known to perform the correction, we chose to perform the sea-level rise forecasts with the uncorrected tidal reconstructions. Additionally, this simplification does not require the same linear stage-discharge relationship to hold as sea-level rises, which would intersect fewer tidal conditions as the harborside water levels increase with sea level. Once the tidal reconstructions for the harborside water levels were made, we added future scenarios of sea-level rise. We used decadal relative sea level forecasts for the nearby Boston Harbor tide gauge (National Oceanic and Atmospheric Administration site 8,443,970) extending from 2020 to 2100 from a joint agency sea level change viewer (https:// geopo rt. usgs. esipf ed. org/ terri aslc/). We selected three global mean sea-level rise scenarios of 0.3, 0.5, and 1.0 m, and we included the low, medium, and high sub-scenarios spanning climate projection uncertainties (Sweet et al. 2017). These sea level forecasts for Boston Harbor were relative to the year 2000 tide gauge water level, so we approximated that the relative increase from 2000 to 2020 for Boston Harbor was consistent with the increase for Wellfleet Harbor and that the discharge from the dike maintained a constant relationship with Wellfleet Harbor water levels in those forecasts. This allowed us to add the sea-level rise projections directly with the tidal reconstructions for the harborside water levels. We used a second-order polynomial interpolation of the decadal sea level projections to a 30-min sampling interval to be consistent with the tidal reconstructions. Together, these steps resulted in nine time series for the harborside water levels with sea-level rise from 2020 to 2100 to set the main boundary condition for the flow model (Fig. 3b).
We ran the hydraulic flow model for each parameter combination with R NP > 0.8 from the 2-year simulations and each sea-level rise scenario with a 30-min time step from 1 January 2020 to 1 January 2100. We developed several metrics to quantify how the water level in the HRE responds to sea-level rise. First, we used the model outputs directly to provide HRE water level, water volume, and water surface area trajectories through 2100 and calculated moving annual arithmetic means to remove the substantial tidal variability in these values. Second, we used the modeled discharge time series to quantify how well the HRE drains over one tidal cycle. We introduce and have defined "drainability" as a new metric to track the draining or impounding behavior of a diked system, but the drainability alone does not quantify the volume stored beyond the drainage capacity of a given tidal cycle. To quantify this drainability, we separated periods of flow into and out of the HRE, representing flood and ebb tide conditions, and we integrated the discharge to yield the volume of water exchanged within each period. The water volumes contributing to the drainability included both flow through the dike structures and freshwater discharge from the contributing watershed (i.e., q in ). Starting with the first flood tide in the time series, we calculated the difference between the flood tide volume entering the HRE with the following ebb tide volume discharged from the HRE. If the ebb tide conditions discharge more water than entered the HRE during the preceding flood tide, we labeled that tidal cycle as drained. We quantify the long-term drainability behavior of the HRE from 2020 to 2100 by calculating the proportion of tides within a year that are drained. Thus, an annual drainability value of one indicates that all tides were net draining, and an annual drainability of zero indicates that, in that year, none of the ebb tides could effectively drain the water flowing into the estuary over the preceding flood tide. Decreasing drainability values track incipient and increasing impounding behavior for a diked system. Annual drainability values of zero establish the complete impoundment of a diked system. We additionally demonstrate the application of the hydraulic flow model by testing how alternative infrastructure affects the response of the HRE to sea-level rise. These tests illustrate how the model could be used to develop water-level based restoration or infrastructure plans. We altered the number of sluices and flap gates in the forecast models, keeping the construction and properties of the existing infrastructure constant. We used a single parameter combination based on the best fit of the sensitivity testing (R NP = 0.97). We ran the models with the infrastructure introduced in 2020, although the infrastructure could have been introduced at any time in the forecasts. We explored the drainage and water storage behavior of the HRE with only a sluice or only flap gates, increasing to five flap gates with one sluice, increasing to three sluices with two flap gates, and adding one sluice and one flap gate to the existing infrastructure.

Parameterized Present-day Model
Sensitivity testing for the HRE dike flow control structures identified the sluice free-flow flow coefficient, C, during the flood stage as the most sensitive model parameter of all parameters tested (Fig. 4a). The upstream flow input into the HRE, q in , was the next most important parameter followed by the sluice free-flow headloss terms (i.e., H loss_max and D HL ). These results show the importance of the two hydrologic inputs to the HRE, q in and C during flood stages, with the parameters controlling the drainage contributing less to the diked water levels. Our sensitivity results are site-dependent and controlled by the flow control structure construction (e.g., culvert dimensions, gate properties, and factors influencing C in Eq. (1)) as well as the hydrologic setting. Next, we confirm that the hydraulic model can accurately reproduce observed water levels within the HRE. A total of 249 30-day sensitivity testing models out of the 38,912 resulted in R NP > 0.9 (Fig. 4b), with 1035 with R NP > 0.8 and 27 with R NP > 0.95. Of those 249 models, 45 models resulted in a R NP > 0.8 for the 2-year, 30-min time step present-day models. Models with R NP < 0.8 for the 2-year simulations generally predicted overly high water levels in the HRE with insufficient drainage capacity during low tides. The 45 higher efficiency model results matched the amplitude and timing of observed HRE water levels with residuals near zero (i.e., modeled values subtracted from observations) (Fig. 5). The median interquartile range (25th to 75th) for the residuals across the 249 high R NP models was − 0.02 m (i.e., overprediction) to 0.11 m (i.e., underprediction). For 166 of the high R NP models, flow instabilities resulting from too much drainage to the harbor led to short-lived (e.g., one tidal cycle) and unrealistically low modeled HRE water levels with residuals reaching > 1.5 m (Figs. 3a and 5a). Such excesses of outflow could be reduced by using a shorter time step or by enforcing a minimum water level allowed for outflow equal to the lowest invert elevation. A subset of 69 of the 249 high R NP models resulted in residuals always < 0.5 m.
The maximum overprediction of HRE water levels in each high R NP models had a median of − 0.35 m. Overall, the hydraulic flow model provides useful information and forecasts estuarine water levels accurately for most tidal conditions (Fig. 5). The ensemble of numerous parameter combinations leading to high-efficiency models contained some short-term artifacts that did not influence the longterm water level forecasts. High-efficiency models did not demonstrate a clustering of hydraulic parameter values and instead spanned nearly the full range of the parameters input to the Sobol sensitivity analysis ("Sobol Sensitivity Testing" section). We used the 45 2-year high efficiency parameter combinations in the following analyses to constrain the long-term sea-level rise and alternative infrastructure scenario HRE water level forecasts. Each sea-level rise scenario consisted of 135 models, and 405 sea-level rise models with the existing infrastructure were run from 2020 to 2100.

Sea-level Rise Scenarios
All sea-level rise scenarios for the HRE with the existing infrastructure indicated a conversion from somewhat drained conditions in 2020 to nearly complete impounded  conditions by 2100 (Fig. 6). Across the parameter combinations, the present-day simulated median drainability was ~0.3, signifying that on an annual basis only 30% of the tides today lead to net water removal from the HRE. The 95th percentile for present-day drainability was ~0.8 and the maximum simulated current drainability was ~0.9. The present-day low drainability results show that the existing infrastructure, designed to create net draining conditions in the HRE, are already under performing for a majority of tides. However, the drainability metric does not account for integrated drainage occurring over or across multiple tides. The simulated net annual storage change in the HRE also suggested the existing infrastructure would be incapable of draining the HRE in the future (Fig. 7). The water volume stored in the HRE was simulated to increase by 3-20% by 2030, 7-77% by 2050, and 11-600% by 2100 across the wide range of sealevel rise scenarios examined relative to the mean 2020 HRE volume. These results support the hypothesis that the HRE will become increasingly impounded with sealevel rise if no changes are made to the flow structures or system hydrology.

Alternative Infrastructure Scenarios
With the impoundment of the HRE likely with the existing hydraulic infrastructure in the dike, we tested several alternative infrastructure solutions to improve drainability in the future. Introducing more flap gates with the same construction as the existing flap gates extended the drained conditions within the HRE by up to several decades. Each extra flap gate accounted for roughly a decade extension of drained condition with the uncertainty in sea-level rise subscenarios leading to also about a decade of uncertainty for major losses in drainability (Fig. 8a-c). Additional sluice gates with the same construction as the existing sluice resulted in more rapid and severe impoundment relative to the existing infrastructure ( Fig. 8d-f). None of the alternative infrastructure scenarios tested provides clear maintenance of the existing drainability conditions at the HRE, although the use of five flap gates with the lowest sea-level rise scenario and climate sub-scenario was effectively stable to 2100 (Fig. 8c). More hydraulic structure scenarios could be tested to optimize the drainability of the HRE, including testing more gates or gates with different designs and flow coefficients.

Model Limitations and Opportunities for Enhancement
The implementation of the hydraulic flow model is computationally efficient, allowing rapid and numerous simulations with relatively few input parameters. Conceptually, the model consisted of a prescribed water level on one side of the hydraulic flow structures and the volumetric flow rates through those structures for a water level difference across the structures. The representation of how the diked water levels responded to input flow was simply a stage-storage relationship for the HRE. The stage-storage relationship removed the hydrologic complexities of flow within the estuary, not accounting for observed spatially varying water levels, flow velocities, salinities, or groundwater feedbacks (WHG Inc. 2012; S. M. Smith and Medeiros 2019; Mullaney et al. 2020). The model output within each time step is simply the volumetric flow rate across the dike flow control structure. This flow rate could be used to set a flux boundary condition in more advanced models that aim to understand spatially varying hydrology within the diked system, as was previously done for the HRE (WHG Inc. 2012).
Similarly, we tested the simulation using constant flow from the contributing Herring River watershed. Despite this simplification, the present-day modeled HRE water levels generally performed well, even without seasonal or event-based changes in streamflow in the Herring River. The hydraulic flow model was built to allow a time-dependent source term, which could also sum the streamflow and direct surface runoff to the waterbody of interest with any other direct water inputs (e.g., groundwater discharge, storm sewer inflows, precipitation). This flow input could also be used to simulate the long-term effects of climate change on precipitation or the water balance more generally. Such a study could investigate the feedback between climate change and sea-level rise on diked water levels.
For the HRE, increasing and more extreme precipitation is predicted for the next century (Massachusetts Environmental Policy Act Office 2018). Thus, q in would be expected to rise over our future scenario simulations, unless groundwater pumping from local municipalities reduced groundwater inputs to the HRE. The additional freshwater influx with climate change would likely accelerate the loss of drainability and create more impoundment than forecasted in our static q in scenarios. Shorter-term variability in q in related to the seasonality of precipitation, evapotranspiration, and pumping would influence the drainability calculations on weekly to monthly timescales. Additional feedbacks could occur with sea-level rise, where rising groundwater levels could lead to more baseflow to the HRE (Masterson and Garabedian 2007). Since the average magnitude of q in is less than ~1/10th of the present-day outflow from the existing infrastructure (WHG Inc. 2012), we anticipate q in dynamics and increases caused by climate change to be less influential on drainability than sea-level rise. Therefore, we chose to focus on the sea-level rise scenarios in the absences of these other hydrologic drivers to understand the long-term implications of the hydraulic infrastructure at the HRE. However, changes in the terrestrial hydrologic system and marine hydrodynamics with climate change could create sufficiently different fluxes and water levels to require incorporation in more detailed studies or for other diked or tidally restricted systems.
For the sea-level rise and alternative infrastructure scenarios, we tested how the HRE water budget would respond to specific forcing or structural changes to the system. Additional studies could investigate how incremental or delayed changes to the infrastructure could stabilize drainage and water levels. For example, additional flap gates could be installed every year to counter the time-dependent drainability losses. Adjustable flow control features could be tested and allow dynamic management of opening sizes or number of outlets. Precipitation or storm surge events could also be developed to simulate shorter-term diked system drainage performance, although more complex hydrodynamic models would be needed for constraining the exterior or harborside water levels. Under diked system restoration plans, this modeling framework could also be used to simulate incremental dike opening effects on surface water levels, although only basic salinity estimates could be calculated from the volume of tidal exchanges.
The model was developed to understand how the HRE would change with sea-level rise, but we generalized the model to allow it to be applied to other diked systems. As an open-source Python library (see Code and Data Availability), the model scripts and functions can be altered and expanded upon for other analyses and to incorporate other Overlapping semi-transparent shaded areas result from the ensemble of 45 high R NP parameter combinations applied to each sea-level rise scenario and climate sub-scenario. Bold lines border each shaded area to highlight the forecast boundaries for a sub-scenario flow structure types (e.g., round flap gates). Our intention was to share the simple model as a useful screening tool that could later justify more complex hydrodynamic modeling or observation network development. The application of the model is site-specific, with hydrologic, topographic, and infrastructure-based inputs to parameterize for each system. We provided one example of how to use Sobol sensitivity testing to quantify the uncertainty of the model parameterization, but alternative sensitivity testing or optimization methods could be used to further understand dike flow control structure systems.

Socio-ecological Implications
Understanding how diked coastal hydrologic features respond to sea-level rise remains a challenging and a coupled problem across ecological and socioeconomic scales. Fundamentally, most diked systems today were the product of human socioeconomic and political choices, where environmental and ecological considerations have only been recognized in recent decades (Roman et al. 1984;Dionne et al. 1998). Such human choices will determine how, when, or if diked or otherwise tidally restricted systems are neglected, expanded, retrofit, managed, or otherwise restored to natural conditions in the next several decades. Coastal infrastructure could have been historically designed under the incorrect expectation of long-term hydrologic stationarity or for shorter operational lifetimes than their eventual tenure. The application of modeling the current and future drainability of such systems could identify priority systems for improvement of ecosystem function and/or flood control. The removal or alteration of existing diked systems could also be designed to slowly change hydrologic conditions, such as enhanced tidal mixing and increased salinity, allowing gradual ecosystem changes. For example, managed hydrologic conditions could aim to restore a salt marsh where hydrologic management actions have resulted in freshwater wetlands as displaced ecotones within the environment. In addition, as is planned for the HRE, ecosystem restoration could include initial elevation enhancement (e.g., sediment augmentation) to overcome subsidence and provide sufficient elevation capital to avoid flooded conditions with the removal of flow control structures (Cape Cod National Seashore and Herring River Restoration Committee 2012; Smith et al. 2020). Thus, maintaining existing or historic drainability is not a management goal for the HRE. Modeling changes in drainability could also be used to transition from gray infrastructure to more natural or nature-based strategies for currently diked systems (Sutton-Grier et al. 2015). New coastal infrastructure that enables a return to more natural hydrology may provide pathways to avoid the physiochemical diked ecosystem conversion (Cadier et al. 2020). In the short term, additional infrastructure to allow more drainage only during low tide could maintain drained freshwater conditions. Over the long term, rising low tide elevation with sea-level rise would eventually require raising the barriers and flow structures or require pumping. Alternatively, restoring a tidal connection to diked lowlands can expand the habitat for saltmarsh ecosystems that sequester carbon, providing both protective function and a suite of other ecosystem services (e.g., Miller et al. 2008;Karberg et al. 2018;Janousek et al. 2021). The stored soil carbon can aid saltmarsh soil elevations to keep pace with sea-level rise and subsidence, offering a potentially longer-term solution than infrastructure upgrades.

Conclusions
We used a simple hydraulic flow model to forecast water levels of the HRE in Wellfleet, MA, to 2100 with expected trajectories of sea-level rise. We found that the HRE may already be incapable of draining the tidal, riverine, and groundwater inputs it receives over each tidal cycle. Simulation ensembles forecast increasing impoundment with sea-level rise over the next several decades. To further understand this conversion from drained to impounded conditions, we developed and demonstrated the application of a novel metric for monitoring and forecasting the drainability of diked systems with sealevel change. The drainability quantifies the net drainage or impoundment state occurring over tidal cycles. In our analysis, the ensemble median present-day annual drainability of the HRE was only 30% and decreased to below 10% by 2060. Fewer low tides per year that allow drainage of diked water leads to increasing water levels and net impoundment. Drainability can also be used in diked systems with one-way flow (e.g., tide gates only) to quantify freshwater or upland drainage performance across tides, as well as in natural or built tidally restricted systems. Similarly, the drainability metric could be calculated over different time spans, including lunar cycles or seasons, to provide shorter-term management and planning insights. The drainability metric that we introduced can be used as a normalized criterion for monitoring, managing, and forecasting the behavior of diked lowland systems. Drainability can aid in water level and volume analyses to develop site-specific thresholds for systems facing impoundment with relative increases in sea level.

Code and Data Availability
The dike hydraulic flow model code, supportive functions, and post-processing algorithms are available at github. com/ kbefus/ dike-flow-hydra ulic. All other datasets are either produced by these scripts or are available from the original citations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.