Mesh Development and Topobathy

This study uses a nesting approach that features three levels of computational domains, each discretized with curvilinear meshes. The level 1 ocean-scale domain (L1) spans part of the north Atlantic Ocean, Gulf of Mexico, and Caribbean Sea. While the L1 mesh structure aligns with that of Zhu and others (2023), modifications have been made to refine the mesh specifically at Shinnecock Bay. The level 2 bay-scale domain (L2) (depicted by the green box in fig. 4) encompasses the entire Shinnecock Bay, measuring approximately 24,488 m alongshore and 10,158 m cross shore. Its south boundary extends approximately 4,400 m into the Atlantic Ocean, with a water depth of around 26 m. The north boundary predominantly lies on land but traverses Cold Spring Pond adjacent to Great Peconic Bay and the Shinnecock Canal connecting Shinnecock Bay and Great Peconic Bay. The level 3 local scale domain (L3) (indicated by the blue box in fig. 4) covers the Shinnecock Indian Nation and the eastern part of Shinnecock Bay. Overall, the mesh sizes gradually decrease from 40 kilometers in deep basins and flat shelves (L1) to around 1 m at the study site (L3). Table 1 lists the mesh size ranges within the three levels of computational domains.

Topographic and bathymetric data are obtained from multiple sources: (1) the northeast Atlantic 3 arc-second digital elevation models (DEMs) with about 90 m horizontal resolutions (NOAA National Geophysical Data Center, 1999); (2) NOAA and Federal Emergency Management Agency (FEMA) bathymetry data with horizontal resolutions ranging from 30 to 2 m (FEMA, 2014); and (3) USGS Coastal National Elevation Database topobathy DEM compiled in 2016 with 1 m resolution based on lidar data for topography and 10 m resolution for bathymetry (Danielson and others, 2016).

Model Coupling, Boundary Conditions, and Meteorological Forcing

Delft3D-FLOW and SWAN are one-way coupled to optimize computational efficiency. The spatial-temporal varying water levels and depth-averaged current fields from the Delft3D-FLOW model feed into SWAN for wave simulations; however, the wave fields from SWAN are not sent back to the FLOW model for flow simulations. While in deep waters, the flow fields are less affected by waves; in the surf zone, they are more subject to the influence of wave actions because of (1) enhanced vertical mixing processes attributed to turbulence generated by white capping and depth-limited wave breaking; (2) increased weight of wave-induced bottom shear stress; and (3) wave setup/setdown induced by the gradient of radiation stress. In the nearshore, Zhu and others (2023) demonstrated that two-way coupling slightly outperforms one-way coupling in water level and wave simulations. Moreover, running L3 Delft3D-FLOW for the whole year is computationally demanding; therefore, we use water levels from L2 Delft3D-FLOW as input to L3 Delft3D-WAVE. Given that the effects of boulders on water levels and depth-averaged current fields are not validated due to the lack of field data, the current fields are not activated in L3-WAVE.

At the seaward boundary of the L1 domain, the Delft3D-FLOW model is forced with the astronomical tide with seven principal constituents determined from a tidal constituent database (Mukai and others, 2002). In SWAN, integral wave parameters (that is, H_s, T_p, and θ_p) from the global and regional scale wave model (WW3 model) are enforced along the seaward boundary to ensure that long wave signals are included in L1 to L3 SWAN models (Zhu and others, 2023). The Delft3D-FLOW and SWAN models in the upper-level domain provide boundary conditions, including hourly outputs of water level, current, and wave spectra, for the lower-level or nested domain.

The north boundary of the L2 domain intersects the Shinnecock Canal, an artificially cut and maintained conveyance channel that allows water to flow from Great Peconic Bay to Shinnecock Bay (fig. 1). A gate that prohibits water from flowing in the opposite direction is operated based on the disparity in water levels between Great Peconic Bay and Shinnecock Bay. Water levels are monitored at the USGS stations in Great Peconic Bay (USGS01304562; USGS, 2024a) and Shinnecock Bay (USGS01304746; USGS, 2024b), as depicted in figure 5 (refer to figure 4 for station locations). Due to the difference in tidal phases between the two bays, high tide and low tide events occur at distinct times on each end of the canal. The peak highs in Shinnecock Bay lead the peak highs in Great Peconic Bay by around 5 hours (see upper panel in fig. 5). In Shinnecock Bay, high and low tide stages correspond with rising and falling tides, respectively, in Great Peconic Bay. Conversely, in Great Peconic Bay, high and low tide stages correspond to falling and rising tides, respectively, in Shinnecock Bay. The gate opens when the water levels in Great Peconic Bay surpass those in Shinnecock Bay and closes when the water levels in the two bays are nearly equal. Representative values for the water level difference are 30 and −5 centimeters for gate opening and closing, respectively (Militello and Kraus, 2001). To simulate the timing and one-way flow through the Shinnecock Canal, a composite water level boundary condition, computed as the maximum water levels at two ends of the canal (blue line in fig. 5), is implemented at the north boundary of the L2-FLOW model.

The U.S. Army Corps of Engineers Engineering Research and Development Center’s Coastal and Hydraulics Laboratory deployed a current meter in the Shinnecock Canal in 1998 to monitor the current speed in the canal (Militello and Kraus, 2001). In this study, we locally increased the bottom friction factor (ƒ₀) from 0.02 to 0.04 to align the modeled maximum current speed with the field measurements. Figure 6 displays the depth-averaged velocity from L2-FLOW without composite water level conditions and with composite water level conditions with ƒ₀ = 0.02 and ƒ₀ = 0.04. The figure shows time intervals during which the current speed approaches near zero. The initiation of these intervals corresponds to times of gate closings, while their conclusion aligns with gate openings.

The meteorological forcing comes from two sources. In the L1 and L2 domains, the spatial-temporal varying wind and pressure fields are obtained from the NOAA North American Mesoscale (NAM) Forecast System, whereas in the L3 domain, spatially uniform and temporally varying wind and pressure fields are obtained from observations at National Data Buoy Center (NDBC) station 44017 (fig. 4). Missing observation data in January 2018 is supplemented by data from a nearby station, 44025 (fig. 4). Figure 7 shows that the two stations recorded similar wind speed and direction, and their wind roses have high resemblance. Wang and others (2022) applied an artificial neural network technique to address missing wind data at NOAA stations. The same technique can be applied to obtain more accurate estimations for absent wind data.

Treatment of Boulders

The length of the boulder line along the western facing shoreline of the Shinnecock Indian Nation is approximately 722 m. Distances between the boulders range from around 0.3 to 5 m (fig. 2). The L3 computational mesh is still too coarse to resolve the boulders, especially to the north and south of the Shinnecock Indian Nation shoreline. In the L3-FLOW model, boulders are represented as porous plates, which are groups of points in the computational grid where resistance to flow is added to the governing equations to simulate the effects of a porous structure. A quadratic friction term as below is included in the momentum equations to parameterize the force on the flow generated by the boulders:

where c_loss is the energy loss coefficient. Porous plates can only be defined at multiples of 45 degree angles with the grid directions. We assume that each cell with boulders features two porous plates, one in the u direction and one in the v direction. Each porous plate cell is assigned a constant c_loss value that partially allows for the transfer of mass and momentum across the porous plate. To the best of authors’ knowledge, the value of c_loss for boulders is currently unavailable. Additionally, the current reduction behind the boulders is not available, presenting challenges in the calibration of c_loss. In response, we ran L3-FLOW with c_loss = 1.0, 10.0, 40.0, and 100.0. The sensitivity of model results to different c_loss values is demonstrated in the results section.

L3-SWAN is carried out both with and without the presence of boulders to evaluate the impact of boulders on wave height. Boulders are represented as thin walls, dams, and vegetation. Thin walls and dams are obstacles that can reduce wave height, reflect waves, and cause diffractions around their ends. A transmission coefficient (K_t), computed as below (Goda and others, 1967), is applied in the propagation terms of the action balance equation:

where We use the temporally spatially varying H_s,i and F to compute K_t, leaving α and β the only calibration parameters. According to Seelig (1979), α = 1.8 and β = 0.1 correspond to vertical thin walls, while α = 2.6 and β = 0.15 correspond to dams with the slope of 1:3/2. Different from thin walls and dams, vegetation contributes solely to wave damping. The energy dissipation rate due to vegetation is included as an energy sink term in the wave action equation. The time-averaged rate of energy dissipation per unit horizontal area over the entire height of vegetation induced by random waves is (Mendez and Losada, 2004): where Boulder parameters from the lidar survey (fig. 3, average height = 0.79 m, average width = 0.5 m, density = 12.85 boulders/m²) were used to compute boulder-induced wave energy dissipation using equation 4. To summarize, we let b_v = 0.5 m, N_v = 12.85 boulders/m², and h_v = 0.79 m. The drag coefficient C_D is the only model calibration parameter. Wave heights before and after the boulders are not available for calibrating C_D, therefore, we set C_D = 1 to demonstrate the vegetation effects on wave heights. Field data is desired for more accurate model calibration.

In Delft3D-FLOW models, the spatial distributions of Manning’s coefficient come from the ADCIRC model used in previous surge studies (Dietrich and others, 2011). The time steps for Delft3D-FLOW models are set as 2, 0.05, and 0.005 minutes for L1, L2, and L3 domains, respectively, to ensure model stability. In SWAN models, the wind input applied in the saturation-based parameterization follows Yan (1987), and the Joint North Sea Wave Project (JONSWAP) coefficient is set as 0.038 square meter per cubic second (m²/s³) (Hasselmann and others, 1973). The depth-induced wave breaking is simulated using a constant breaker index of 0.73. The time steps of SWAN models in L1, L2, and L3 domains are all set as 1 hour.

Model Performance Metrics

To assess the model performance, the following statistical parameters are utilized: normalized root-mean-square difference (NRMSD), coefficient of determination (R²), and scatter index (SI). The definitions of these parameters are provided below:

The root-mean-square difference (RMSD) is defined as

where The RMSD is normalized by the difference between the maximum and minimum observations to compute NRMSD.

The coefficient of determination R² is defined as

where the sum of squares of residual RSS is computed as $RSS={\displaystyle \sum }_{i=0}^N{\left(O_i-M_i\right)}^2$. Ō is the average of the observations. The total sum of squares TSS is computed as $TSS={\displaystyle \sum }_{i=0}^N{\left(O_i-\overline{O}\right)}^2$. R² = 1 indicates perfect agreement between the modeled and measured variables. R² = 0 indicates that the model results are as good as random guesses around the mean of the observations. R² is different from the square of Pearson correlation coefficient r, which measures the linear correlation of the modeled and measured variables.

The scatter index (SI) is defined as

Model Validation

The L1-SWAN model is validated against measurements from the NDBC station 44017, which is approximately 39,650 m southeast of Shinnecock Inlet (see figure 4 for the buoy location). The modeled wave parameters, including H_s, T_p, and θ_p, agree well with measurements (figs. 8, 9). The model performance metrics are summarized in table 2. Despite the relatively coarse grid at the NDBC station 44017 location, which is nearly 9 kilometers, the measured and modeled probability histograms of T_p align well (fig. 10). The probability histograms of T_p indicate a peak around 6 seconds. The presence of longer period swell waves with T_p greater than 12 seconds or shorter period waves with T_p less than 2 seconds is much less likely. The most likely T_p falls within the range of 4–9 seconds, characteristic of typical wind-driven seas. The absence of wave measurements within Shinnecock Bay prevents model validation of local-scale wave parameters. While we are unable to directly validate the L2 and L3 wave models, the strong agreement in wave parameters at NDBC station 44017 from L1-SWAN demonstrates that accurate wave boundary conditions can be transferred from L1-SWAN to L2-SWAN. The complex nearshore wave processes, such as wave breaking and wave shoaling, are considered in both L2-SWAN and L3-SWAN models along the shoreline of Shinnecock Bay. The local-scale wave dynamics and parameters in Shinnecock Bay are simulated reasonably by the Delft3D-FLOW and SWAN modeling system. Zhu and others (2023) used similar nested domains and coupled Delft3D-FLOW and SWAN models to simulate wave dynamics in upper Delaware Bay with oyster reef-based living shorelines. Local-scale field wave measurements along the Shinnecock Indian Nation shoreline could further validate the wave model and confirm the validity of this modeling system in predicting local-scale wave dynamics along the Shinnecock shoreline.

The modeled water level (η) from L2-FLOW is validated against measurements from the USGS station 01304746, which is close to the Ponquogue Bridge in the middle of Shinnecock Bay (see fig. 4 for station location). The modeled water level aligns well with measurements (fig. 11). The model metrics are summarized in table 2. The good agreement in water levels of L2-FLOW at USGS station 01304746 substantiates that accurate water level boundary conditions are fed into L3-FLOW from L2-FLOW, although there is no water level data near the Shinnecock Indian Nation shoreline for the L3-FLOW model calibration and validation. It should be noted that field measurements of current velocity along the shoreline in Shinnecock Indian Nation are not available; therefore, comparisons between simulated and observed current velocities for further model validation is impossible.

Wave Climate From L2-SWAN

Figure 12 shows the elevation map of the eastern part of Shinnecock Bay, where the offshore elevation maintains a nearly constant level of around −3.5 m (North American Vertical Datum of 1988). We selected four numerical output stations (A–D in fig. 12) along the −3.3 m (North American Vertical Datum of 1988) elevation contour and plotted the corresponding wave roses. Figure 12 depicts that the majority of offshore waves approach the Shinnecock Indian Nation shorelines at normal wave angles. The yearly average wave height in 2018 is around 0.19 m, with a standard deviation of around 0.14 m. The maximum modeled wave height in 2018 reaches 0.65 m. The wave power (P) is computed as below:

where P_x and P_y are the energy transport components computed internally in SWAN based on the variance density spectra. The yearly averaged P at output stations A–D are 49, 50, 48, and 44 watts per meter (W/m), respectively.

Boulder Effects on Currents

To illustrate the effects of boulders on currents, a baseline run without boulders is carried out in L3-FLOW. Numerical output stations are placed around the boulders as shown in figure 2. Figure 13 shows the magnitude of the daily averaged depth-averaged current (|U|) at upwave and downwave output stations corresponding to different c_loss. In general, the downwave |U| decreases as c_loss increases. Following a linear regression between the upwave and downwave daily averaged currents (U), we compute the reduction rate of the daily averaged currents as 1 − m, where m is the slope of linear regression. The sensitivity of the reduction rate to c_loss is illustrated in figure 14. At the output stations near the tidal ponds, a significant reduction in U is observed in the baseline run, and further reduction in U is induced by the presence of boulders. In Delft3D-FLOW, the vertical position of the porous plate is specified in terms of the model layers that it occupies. In this study, only one vertical layer is used; therefore, the boulders will always result in momentum loss regardless of the submergence. Additionally, the calibration of c_loss relies on the measured data in water level and depth-averaged velocity around the boulders. Due to the lack of field measurements, this study only provides a modeling analysis of the effects of boulders on currents. Simulation results show that the current velocity at the output stations could reach as high as 0.18 meter per second (m/s) without boulders as shoreline protection structures but can be reduced to less than 0.02 m/s (fig. 13). The current velocity reduction rates when c_loss = 100 at the case with boulders compared to the case without boulders increase from less than 20 to larger than 70 percent at the P3 and P4 stations, from less than 30 to larger than 80 percent at the P5 and P6 stations, from less than 60 to larger than 80 percent at the P7 and P8 stations, and from less than 30 to larger than 85 percent at the P9 and P10 stations, respectively (fig. 14).

Boulder Effects on Waves

Model results indicate variability in longshore wave heights. H_s is generally higher near the north end of the boulders (P3 and P4 stations; refer to fig. 2 for output station locations) and lower near the tidal ponds south of the cemetery (P7 and P8 stations). Figure 15 illustrates wave height changes at the upwave and downwave sides of the boulders. Generally, wave heights decrease at the downwave side of the boulders. Wave height at the output stations in the upwave side of the boulders range between 0.35 and 0.5 m, whereas at the downwave side it can be reduced to less than 0.1 m (fig. 15). Both the dam and thin wall treatments of boulders in SWAN show similar H_s variations. The average reduction rates of H_s are 52, 15, 63, and 15 percent at P4, P6, P8, and P10 stations, respectively (fig. 15). The vegetation treatment of boulders results in generally smaller H_s reduction. The reduction rates are 16, 13, 63, and 26 percent at P4, P6, P8, and P10 stations, respectively (fig. 15). In addition to wave-height reduction, wave amplification is also observed in the model results. To illustrate the effects of boulders on wave transformation, a baseline run without boulders is carried out in L3-SWAN. The model results considering boulders are then compared with baseline model results (see fig. 16). At the north end, the dam and thin wall treatments of boulders result in a 44 percent reduction in H_s, while the vegetation treatment (treating boulders as vegetation) of boulders only leads to a 2 percent reduction in H_s. At the south end, the dam, thin wall, and vegetation treatments of boulders lead to 5, 2, and 2 percent reductions in H_s, respectively. Near the two tidal ponds, boulders do not cause H_s changes on average.

The yearly averaged wave power at the upwave side of the boulders varies from 2 (P7 station) to 35 (P3 station) W/m, with the smallest yearly averaged wave power found near the tidal pond south of the cemetery. In the baseline run, the yearly averaged wave power is reduced to 25, 3, 0.2, and 10 W/m at P4, P6, P8, and P10 stations, respectively (fig. 17). Treating boulders as dams or thin walls leads to additional wave power reductions of 13, 0.02, 0.01, and 1 W/m at P4, P6, P8, and P10 stations, respectively (fig. 17). Treating boulders as vegetation leads to additional wave power reductions of 3, 0.2, 0.01, and 4 W/m at P4, P6, P8, and P10 stations, respectively (fig. 17). The annual wave power reduction around the boulders is consistent with the wave height variations, with boulders causing significant annual wave power reduction near the north end, nearly unchanged annual wave power near the two tidal ponds, and small annual wave power reduction near the south end of the boulders.

Effect of Living Shoreline Structures on Wave Dynamics

Simulation results show that wave height decreases at the downwave side of the boulders along the Shinnecock Indian Nation shoreline, indicating that the boulders are generally effective in attenuating wave energy. There are large local variabilities in wave height and wave power reduction from upwave output stations to downwave output stations. Reductions in wave height tend to be larger at the northern end of the boulders where larger wave heights are found when compared to the southern end of the boulders. The variations in wave height reduction result from a combined effect of bathymetric changes around boulders. This is consistent with findings from literature on oyster reef-based living shoreline structures that wave attenuation by living shoreline structures is dependent on water depth, wind speed, wind direction, and local bathymetry (Wiberg and others, 2019; Zhu and others, 2020; Wang and others, 2021, 2023).

Model results also indicate wave amplification—this is consistent with the field measurements associated with oyster reef-based living shoreline structures in Gandys Beach (Zhu and others, 2020). Zhu and others (2020) revealed that constructed oyster reefs do not necessarily induce wave-height reduction, and wave amplification occurs due to wave shoaling, wave diffraction, and wave focusing. It is found that treating boulders as dam, thin wall, or vegetation leads to different magnitudes in reductions of wave height and wave power, demonstrating that different shoreline structures have various effects in wave attenuation. The vegetation treatment of boulders results in smaller wave-height reduction because it does not induce wave breaking when the boulder freeboard approaches zero.

Effect of Living Shoreline Structures on Current Dynamics

Simulation results show that current velocities from bayside of the boulders to leeside of the boulders along the Shinnecock Indian Nation shoreline are reduced more with the boulders compared with the case without the placement of boulders. Moreover, reductions in current velocities along the Shinnecock Indian Nation shoreline tend to vary with locations on the shoreline. Compared with the no-boulder case, current velocity reductions with the boulders are found to be larger at the output stations that are located on the sections far away from the large boulder gaps (for example, P3 and P4 stations and P9 and P10 stations) than those at the stations at or close to the boulder gaps (for example, P5 and P6 stations and P7 and P8 stations) (fig. 2). These results indicate that the boulders are effective in changing tidal circulation and reducing current velocity, and that there are local variabilities in a boulder’s capacity in current velocity reduction. This finding is consistent with previous living shoreline velocity studies (Wang and others, 2021; Salatin and others, 2022; Wang and others, 2023). Current velocities are found to decline at leeside of the constructed oyster reefs along the Gandys Beach (Wang and others, 2021; Salatin and others, 2022) and Chincoteague (Wang and others, 2023) shorelines, and such velocity reductions by the oyster reefs vary in space and time. Current velocity reductions along the living shorelines are complex and dependent on wind and wave conditions, tidal phase (ebb versus flood), location of measurements relative to the structures, and geomorphological features (Wang and others, 2023).

In this study, although there are no field data to determine which c_loss is representative of the field condition, current velocities in the leeside of the boulders along the Shinnecock Indian Nation shorelines are found to be reduced compared with the bayside of the boulders at the output stations. This result is different from previous studies on current velocity change along the living shorelines using oyster reef-based structure. For example, in Gandys Beach, Virginia (Salatin and others, 2022), and Little Toms Cove, Virginia (Wang and others, 2023), velocities behind the oyster castles can be increased at locations near the boundary between the castles and open water, especially during ebb tides. The discrepancies in structure-induced velocity changes between this study and previous studies can be attributed to the difference in structure configuration as well as the absence of wave-driven currents in this modeling study. The circulation pattern may be altered if the gaps between the oyster castles are larger than the gaps between the boulders (Wang and others, 2021).

Implications for Coastal Management and Restoration Projects

This study demonstrates the offshore wave climate near Shinnecock Indian Nation and quantifies the effects of the existing large boulders on wave attenuation through process-based numerical models. Along the Shinnecock Indian Nation shoreline, wave-height reduction varies along the boulder-protected shoreline, indicating that some sections of shoreline have not benefitted from the boulders. Adjusting the configuration of the existing boulders—for example, by increasing the crest height or the density—at these shoreline sections could help to reduce erosion. Another method for structure adjustment is the utilization of oyster reef-based materials such as constructed oyster castles, oyster shell bags, and cages. This is based on the improvement of current velocity by the construction of the boulders along the Shinnecock Indian Nation shoreline, which is one critical habitat suitability metric for optimizing oyster larval transport, recruitment, and oyster growth. A current velocity threshold of 0.15 m/s is used for optimal oyster growth in oyster reef-based living shoreline assessment (for example, Salatin and others, 2022) and lower thresholds of 0.07 and 0.1 m/s are also determined from oyster growth experiments (for example, literature review in Wang and others [2017b]). In this study, simulated current velocities are found to decrease from larger than 0.15 m/s in the bayside of the boulders (upwave output stations) to less than 0.1 m/s in the leeside of the boulders (downwave output stations).

Other habitat suitability metrics, such as salinity and temperature, could also be considered when applying oyster reef-based structures (with boulders) for enhancing the effectiveness of the living shoreline structures. Seagrasses planted around/near the boulders can also be used to make boulders more effective as living shoreline structures. A numerical modeling study for living shorelines in the Chesapeake Bay using Delft3D-FLOW and SWAN models found that erosion can be drastically reduced when combining riprap, salt marsh, and submerged aquatic vegetation as the living shoreline structures (Vona and others, 2021). The Shinnecock Indian Nation also considered other living components including the construction of oyster reefs built from shell or concrete and the planting of seagrass (for example, eelgrass) to enhance shoreline protection and restoration (for example, Anchor QEA, LLC, The Nature Conservancy, and Fine Arts and Sciences, LLC, 2019). Directly south of the Cuffee Beach parking area, old oyster cages have been placed to attenuate wave energy from large coastal storms, reduce erosion, and facilitate sand accretion. Oysters have been farmed on the eastern side of Shinnecock Bay by the Shinnecock Indian Nation for many years and were introduced in the western part of the bay as part of a research and restoration program in 2017 (Ahn and Ronan, 2020). Regardless of the types of living shoreline structures that may be used along the Shinnecock Indian Nation shoreline, the model framework developed in this study can be leveraged for optimal design of nature-based solutions, which can guide where to install living shoreline structures and how to determine their optimal size.

Data and Knowledge Gaps

The reduction in wave height and wave power along the Shinnecock Indian Nation shoreline on a local scale is determined by the L3-SWAN model. Due to the lack of measured wave and velocity data around the boulders, L3-SWAN model simulations could not be validated. This study provides a modeling analysis of the potential effects of boulders or similar structures on currents and waves. The lack of hydrodynamic and wave data hinders the calibration and validation of wave and flow models, limiting our ability to assess and predict local-scale wave dynamics accurately. Addressing this data gap could help to enhance the reliability and effectiveness of hydrodynamic and wave models and could provide valuable insight into the local wave-climate and estuarine circulation—these validated models can help to facilitate informed decision making by the Shinnecock Indian Nation and various stakeholders in the region.

This numerical study reveals that wave and flow measurements near the Shinnecock Indian Nation shoreline, particularly at the coastal restoration project site (for example, boulders and beaches), can help to enhance the understanding of wave climate and coastal processes and calibrate and validate wave models (for example, SWAN in this study). In addition to analyzing wave-height variations from offshore to nearshore, and around the structures, wave-spectral analysis could be conducted to investigate wave-spectral variations. A similar analysis was completed for the oyster reef-based living shoreline near Gandys Beach (upper Delaware Bay) by Zhu and others (2020), which indicated that the energy transfer from the primary waves to the high harmonics after waves propagates over the submerged oyster reefs and revealed that swell energy originated from the Atlantic Ocean can penetrate oyster reefs without any dampening. Current measurements are lacking to calibrate and validate circulation models (for example, Delft3D-FLOW in this study). Boulder dimensions, including the cross-sectional areas, heights, gap width, and population density, are lacking for boulder implementations in numerical models. The nearshore bathymetry and topography are outdated, and thus, surveys and updates in numerical models could help to achieve higher model accuracy. The field data can also be employed in data-driven models, such as the bagged regression tree algorithm (Wang and others, 2022; Zhu and others, 2023) to improve the accuracy of physics-based models at locations where field measurements are available.

Further Research Opportunities

This study found local variabilities in current velocity and wave height and wave energy along the Shinnecock Indian Nation shoreline, and therefore, field measurements of current velocity and wave parameters at multiple locations in both the bayside and leeside of the boulders over the storm periods could help improve the accuracy of the hydrodynamic and wave modeling. Since the main goal of living shorelines is to protect shoreline from wave and current induced sediment erosion and shoreline retreat, the in situ measurements of sediment movement including sediment transport, erosion, deposition, and trapping by vegetation from offshore to onshore with and without living shoreline structures (large boulders for Shinnecock Indian Nation shoreline in this study) could help improve the understanding of sediment dynamics. With field sediment measurements, sediment budget can be determined to inform effective management of sediment along the Shinnecock Indian Nation shoreline and beaches. By integrating sediment dynamics with current and wave dynamics along the Shinnecock Indian Nation shoreline, direct relationships between wave energy and shoreline retreat rates can be explored (Marani and others, 2011; Wang and others, 2023; Zhu and others, 2023). With field sediment data, models of geomorphology for shorelines, beaches, and salt marshes in Shinnecock Indian Nation can be developed to assess and predict the morphological changes in Shinnecock Indian Nation tidal wetlands including shorelines under future storms and sea level rise (Wang and others, 2017a). The modeling system developed in this study can serve as a tool for simulating and analyzing the impacts of major tropical storms and hurricanes. Future studies may utilize this system to better understand the dynamics of extreme weather events and the effectiveness of boulders for shoreline protection and enhance preparedness and response strategies.

Improvement of water quality and habitat for fish and wildlife species and other ecological services are the ecological goals of living shorelines. Oyster reef-based structures and seagrass are considered in Shinnecock Indian Nation living shoreline protection and restoration (Shinnecock Indian Nation, 2013). Field measurements of water quality parameters such as salinity, temperature, dissolved oxygen, total suspended solids, and chlorophyll-a concentrations can help to determine if the habitat suitability for oysters and seagrass species in the estuaries could be improved by implementation of the living shoreline structures (Wang and others, 2023). Field water quality measurements and process-driven water quality models for the Shinnecock Indian Nation estuary can be coupled with the hydrodynamic and wave models (Delft3D-FLOW and SWAN in this study) and validated and used to predict habitat changes in the Shinnecock Indian Nation estuary with and without living shoreline structures under future storms and sea level rise conditions.