Sediment Exchange Between the Created Saltmarshes of Living Shorelines and Adjacent Submersed Aquatic Vegetation in the Chesapeake Bay

Rising sea levels and the increased frequency of extreme events put coastal communities at serious risk. In response, shoreline armoring for stabilization has been widespread. However, this solution does not take the ecological aspects of the coasts into account. The “living shoreline” technique includes coastal ecology by incorporating natural habitat features, such as saltmarshes, into shoreline stabilization. However, the impacts of living shorelines on adjacent benthic communities, such as submersed aquatic vegetation (SAV), are not yet clear. In particular, while both marshes and SAV trap the sediment necessary for their resilience to environmental change, the synergies between the communities are not well-understood. To help quantify the ecological and protective (shoreline stabilization) aspects of living shorelines, we presented modeling results using the Delft3D-SWAN system on sediment transport between the created saltmarshes of the living shorelines and adjacent SAV in a subestuary of Chesapeake Bay. We used a double numerical approach to primarily validate deposition measurements made in the field and to further quantify the sediment balance between the two vegetation communities using an idealized model. This model used the same numerical domain with different wave heights, periods, and basin slopes and includes the presence of rip-rap, which is often used together with marsh plantings in living shorelines, to look at the influences of artificial structures on the sediment exchange between the plant communities. The results of this study indicated lower shear stress, lower erosion rates, and higher deposition rates within the SAV bed compared with the scenario with the marsh only, which helped stabilize bottom sediments by making the sediment balance positive in case of moderate wave climate (deposition within the two vegetations higher than the sediment loss). The presence of rip-rap resulted in a positive sediment balance, especially in the case of extreme events, where sediment balance was magnified. Overall, this study concluded that SAV helps stabilize bed level and shoreline, and rip-rap works better with extreme conditions, demonstrating how the right combination of natural and built solutions can work well in terms of ecology and coastal protection.


INTRODUCTION
Human populations are concentrated in the coastal zone, with three-quarters of the global population living within 50 km of the sea and 50% of the US population living within 50 miles of the sea. More than one-third of the US gross national product is generated in the coastal zone (Marra et al., 2007). Recent catastrophic events such as Hurricanes Katrina in 2005, Sandy in 2012, and Florence in 2018 have shown that coastal communities are at great risk of coastal inundation caused by storm surges and sea-level rise (Li et al., 2020).
Traditional coastal engineering interventions such as breakwaters and sea walls are increasingly challenged by these environmental changes, as their preservation may become unmaintainable (Temmerman et al., 2013). Moreover, as artificial structures, they do not take ecological aspects into account and result in generally, but not exclusively, negative ecosystem impacts (Bilkovic and Mitchell, 2013). In contrast to these interventions, the "living shoreline" technique takes ecology into account by incorporating natural habitat features, such as saltmarshes, into shoreline stabilization. The approach aims to provide the same erosion-control functions of armored structures, while also maintaining the ecological benefits of nature-based solutions (Davis et al., 2015;Scyphers et al., 2015;Gittman et al., 2016a). In this study, we defined a "living shoreline" as a narrow marsh fringe with or without adjacent structures (Burke et al., 2005). Recent studies of living shorelines have highlighted the importance and effectiveness of these nature-based solutions in providing ecosystem services and enhancing coastal resilience by reducing wave energy and facilitating sedimentation (Currin et al., 2010;Manis et al., 2015;Sutton-Grier et al., 2015;Palinkas and Lorie, 2018;Bolton, 2020;Safak et al., 2020).
Saltmarshes and submersed aquatic vegetation (SAV) are key components of the ecosystems in intertidal and subtidal regions, respectively. Saltmarshes are well-recognized for supporting critical ecosystem functions and services such as biological production, nutrient cycling, water quality, and nursery habitats (for examples, see Boorman, 1999, Barbier et al., 2011Grabowski et al., 2012). Saltmarshes also help protect shorelines against erosion, since they can work as natural barriers that dissipate wave energy and reduce sediment re-suspension. They also promote sediment deposition through trapping and in situ production, which provides a mechanism by which the marsh platform can gain elevation to sustain itself in the face of sea-level rise (SLR) (Bricker-Urso et al., 1989;Schuerch et al., 2018). On the other hand, SAV provides food, shelter, and breeding areas for waterfowl, fish, invertebrates, and many other species of aquatic life (Catling et al., 1994;Heck et al., 1995;Noordhuis et al., 2002). It is also widely recognized as an important habitat and indicator of water quality in large rivers and estuaries (Barko et al., 1991;Carter et al., 1994;Short and Wyllie-Echeverria, 1996). Submersed aquatic vegetation also has a well-established impact on water flow and sediment dynamics, wherein flow velocity is significantly reduced within plant stands that facilitates sediment deposition Wharton et al., 2006).
Although the beneficial effect of living shorelines on supporting biodiversity and related ecosystem services is welldocumented in the literature (for example, see Gittman et al., 2016b), the effect on adjacent benthic communities, such as SAV, is still unclear. In particular, the synergies between the created saltmarshes of living shorelines and SAV are not well-understood. Both plant communities trap sediment, which is critical for their resilience to environmental change. However, how sediment might be exchanged between the two communities has rarely been investigated (Donatelli et al., 2018;Zhu et al., 2021). A sediment budget is a useful tool for examining this exchange, as it accounts for both the losses and gains of sediment and can help assess whether a planned restoration of a tidal marsh has a high probability of success (Ganju et al., 2013(Ganju et al., , 2019. In this study, we examined the interaction between the created saltmarshes of living shorelines and SAV in adjacent waters in terms of sediment transport. We used a numerical modeling approach with an aim to provide a sediment budget that encompasses both intertidal and subtidal regions, building on previous work to examine the effect of artificial structures such as breakwaters on sediment transport within marsh platforms (Vona et al., 2020).
Beyond the definition of the sediment budget, this work also analyzed different living shoreline configurations to optimize their realization and make them more efficient. A living shoreline composed solely of saltmarshes is more sensitive to the action of external forcing, and its efficiency may be compromised by extreme conditions. The presence of an adjacent SAV bed may help improve the living shoreline efficacy, and the further addition of an artificial structure such as rip-rap could also counter the threat of extreme events.

Model Description
In this study, we used a double numerical approach to investigate the sediment balance between SAV and saltmarshes. First, we developed a model with Delft3D (Deltares, Netherlands) to validate the deposition rate at three study sites in Broad Creek, within the Choptank River (Maryland, USA). After that, through an idealized model, by coupling Delft3D and SWAN (Simulating Waves Nearshore, Delft University of Technology, Netherlands), we studied the interchange of sediments between the two vegetation communities under different wave height, wave period, and basin slope condition.
Delft3D (Roelvink and Van Banning, 1995;Lesser et al., 2004) is an open-source numerical model system that allows the simulation of hydrodynamic flows, wave generation and propagation, sediment transport, and morphological changes. Its main advantage is that its hydrodynamic and morphodynamic modules are fully coupled, so when the bed topography changes, the flow field adjusts in real-time. The FLOW module (Delft3D) then performs hydrodynamic, sediment transport, and morphological changes by discretizing the related equations on a 3D curvilinear finite-difference grid, solved by an alternating direction implicit scheme. For our simulations, we used the 3D formulation of the hydrodynamic and morphodynamic models implemented in Delft3D. The generation and propagation of random and short-crested waves in shallow and deep waters were computed by SWAN, by taking into account processes such as wave-wave interactions, wave refraction and dissipation (which included bottom friction; Hasselmann et al., 1973), and wave breaking (Battjes and Janssen, 1978). Below, we present the essential governing equations for the model. Further details can be found in the study by Lesser et al. (2004). For notations, refer to Table 1.

Governing Equations
Delft3D solved the mass-balance equation in Cartesian coordinates for an incompressible fluid, under the assumption of shallow water and the Boussinesq approximation: where U, V, and W are the averaged fluid velocity along x, y, and z directions, respectively. The horizontal and vertical momentum equations for unsteady, incompressible, and turbulent flows are: where p is the fluid pressure, f is the Coriolis parameter, ρ is the density, τ xx are fluid shear stresses, and g is the gravity acceleration. Because of the shallow-water assumption, the vertical momentum equation was reduced to the hydrostatic pressure equation. The vertical eddy viscosity was computed by the standard k-ε closure model (Rodi and Scheuerer, 1984). The suspended sediment transport is calculated by solving the three-dimensional advection-diffusion equation: where C is the mass concentration of the sediment fraction, ε s is the eddy diffusivities of the sediment fraction, T s is an adaptation timescale, C eq is the local equilibrium depth-averaged suspendedsediment concentration. In case of erosion and deposition, the exchange between the water column and the bed for cohesive sediments was calculated with the Partheniades-Krone formulation (Partheniades, 1965), while the Van Rijn method was used for the non-cohesive fraction (Van Rijn, xbib1993). As already mentioned, the evolution of wave motion was computed by SWAN by solving the spectral action balance equation: Frontiers in Marine Science | www.frontiersin.org where the left-hand side is the kinematic part of the equation. N represents the action density. The first term represents the local density variation over time. The second and the third ones represent the action propagation in geographical space (along the x-y direction with velocity c x and c y , respectively). The fourth term describes the relative frequency shifting due to variations in depths and currents, while the fifth is the depth-induced and current-induced refraction. Variables c σ and c θ are the velocity propagations in the spectral space. The S in the right-hand side of the equation is the non-conservative source/sink term that represents physical processes such as generation, dissipation, or wave energy redistribution. σ is the wave frequency.

Vegetation Model
The impact of vegetation on hydrodynamics was modeled as an effect on bed roughness and flow resistance. In Delft3D, this was accomplished by using the formulation in the study by Baptist (Baptist, 2005), which modeled vegetation as rigid cylinders characterized by stem diameter (D), height (H v ), drag coefficient (C D ), and density (m). The expression of the Chézy coefficient has been derived in the study by Baptist et al. (2007), building on the results of a 1DV k-ε turbulence model developed in the study by Uittenbogaard (2003), which solved a simplification of the 3D Navier-Stokes equation for horizontal flow conditions with additional assumptions to include the effect of vegetation in the k-ε turbulence closure. The vertical flow velocity profile was assumed to be divided into two zones due to the presence of vegetation: constant (u v ) inside the vegetated patch and logarithmic above.
In the case of fully submerged vegetation, the total shear stress τ t was given by the sum of the bed shear stress τ b and the component due to the vegetation τ v : where ρ is the water density, g is the gravitational acceleration, h is the water depth, i is the water surface slope, C b is the drag coefficient of the vegetation, u v is the uniform velocity component, C D is the bottom Chézy coefficient, H v is the vegetation height, n is the number of stems per unit area, D is the stem diameter, and u v is the uniform velocity component defined as: The vegetated bed bottom shear stress was given by: defined as a function of the total shear stress and the reduction factor f s , which is obtained by replacing Equations (11) in (8): By combining Equations (7) and (12), the vegetated bed bottom shear stress became: where the Chézy friction C rs value is given by: in which k is the Von Karman constant (k = 0.41). At the transition from submerged to emergent vegetation, the second term of the equation becomes zero.
Baptist's equation has been largely evaluated through field data and laboratory experiments and the predicted results have been compared by several studies with experimental data, finding a good fit (Arboleda et al., 2010;Crosato and Saleh, 2011).

Model Set Up
The first numerical approach consisted of modeling the Broad Creek area within the Choptank River. The sediment deposition rates measured in three study areas within the creek, namely, Oppenheim (OPP), Hatton Garden (HG), and Ruesch (RU) (Figure 1) by Bolton (2020) were used as validation data for the model. The numerical domain was composed of 266 cells in the x-direction and 565 in the y-direction, each ∼15 m × 15 m. Bathymetry was obtained from the National Ocean Service 30 m resolution Digital Elevation Model (https://ngdc.noaa.gov/ mgg/bathymetry/estuarine/index.html) of the National Oceanic and Atmospheric Administration and was manually adjusted at the validation sites by adding vegetation to replicate the created saltmarshes of the living shorelines (Figure 2). The initial condition set in the models was a fixed water level at 0.5 m with 5 m of an erodible bed level, composed of 50% sandy and 50% muddy sediment. Non-cohesive sediments were characterized by a specific density of 2,650 kg/m 3 , a dry bed density of 1,600 kg/m 3 , and a mean diameter of 100 µm, while characteristics of the cohesive sediment were chosen in agreement with values provided in the study by Berlamont et al. (1993), wherein specific density was 2,650 kg/m 3 , dry bed density was 500 kg/m 3 , and setting velocity was 0.25 mm/s. In the vertical direction, one single layer was used. The initial sediment concentration was equal to 0.03 kg/m 3 , as found in the study by Ensign et al. (2014) in the oligohaline part of the Choptank River. On the south, west, and east boundaries, water-level variations and suspendedsediment concentrations were inserted. The water level was taken from the Chesapeake Bay Operational Forecast System (CBOFS) at the Choptank River entrance station (just outside the entrance to Broad Creek), from October 18 to November 3 for the years 2017 and 2018. We used two different inputs of water levels to refer to the same timeframe of sediment core samples used in the study by Bolton (2020). The sediment concentration was set at 0.06 kg/m 3 (Ensign et al., 2014). On the north1 and north2 borders, Neumann conditions (normal water level derivative equal zero) were inserted. The typical wind conditions within Broad Creek came from S-S/W with maximum intensity at ∼10 m/s (Windfinder https://it.windfinder.com). The model was therefore forced by the uniform wind coming from the south with an intensity of 20 m/s (dissipative model effects reduce wind intensity around 10 m/s) in order to propagate the sediments into the creek. The input water level was validated with the water level extracted from the CBOFS at a point near the study sites (Figure 2).
The bottom stress was modeled with the formulation in the study by Chézy using a constant value equal to 55 m 1/2 /s. The suspended-sediment eddy diffusivities were a function of the fluid eddy diffusivities and were calculated using a horizontal large eddy simulation and grain settling velocity. The horizontal eddy diffusivity coefficient was defined as a combination of the subgrid-scale horizontal eddy viscosity, computed from a horizontal large eddy simulation, and the background horizontal viscosity, which was set equal to 0.001 m 2 /s 2 in this study (Edmonds and Slingerland, 2010;Nardin et al., 2016). To satisfy the numerical stability criteria of Courant-Friedrichs-Lewy, we used a time step t = 6 s (Lesser et al., 2004). Each run simulates the morphological changes of 150 days under constant conditions during the simulation. To decrease the simulation time, we used a morphological scale factor equal to 10, a user device to multiply the deposition and erosion rates in each t. The morphological factor value was chosen to avoid numerical instability. We then multiplied the obtained deposition rate by 2.4 days to calculate the annual deposition.
Starting from the same numerical parameters used in the first modeling approach, we proceeded with the second idealized numerical model, where we varied wave height, wave period, and basin slope. A double nesting grid was used to propagate the wave motion well (Supplementary Figure 1). The outer and coarser wave grid was composed of 51 and 91 cells in the x and y directions, respectively, with a constant resolution of 100 m × 100 m. The flow domain (nested to the wave grid) was 3 km × 2 km, the computational grid was composed of 102 cells in the x-direction and 126 cells in the y-direction, and it was gradually refined from the eastern side (40 m × 40 m) to the western (10 m × 10 m; Figure 3). Sensitivity analyses showed that the resolution of the domain did not considerably impact wave propagation (Supplementary Figure 2). We simulated three different scenarios, which were (1) marsh only, (2) marsh and SAV, and (3) marsh, SAV, and rip-rap.
To incorporate rip-rap into the simulations, we inserted a non-erodible bottom (thickness = 0) along the entire seaward boundary of the marsh to represent the rocks comprising the rip-rap. The initial condition set in the models was a fixed water level at 0.4 m with 5 m of an erodible bed level, composed of 50% sandy and 50% of muddy sediment with the same characteristics described above. The boundary conditions we imposed were the Neumann type for the northern and southern boundaries and a combination of water-level variation, incoming waves, and incoming suspended-sediment concentration for the eastern boundary. In the vertical direction, five non-homogeneous sigma layers were used, decreasing the layer thickness (%) of the local water depth for each layer going down to the bottom.
We varied the wave height (0.2, 0.3, 0.4, 0.5, 0.7, and 1 m), wave period (2, 3, 4, and 5 s), and basin slope (0.8 and 1%, as found by Koskelo et al., 2018 within the Choptank watershed and  at the Virginia Coastal Reserve), while the suspended-sediment concentration was fixed to 0.1 kg/m 3 . Wave parameters (H s and T p ) were selected to simulate waves generated in the Chesapeake Bay (Lin et al., 2002)  and were directed orthogonally to the shoreline. We imposed these values at the eastern boundary of the wave grid. Wave reflection was not accounted for in the wave model, so wave energy was dissipated at the coastline. To satisfy the numerical stability criteria of Courant-Friedrichs-Lewy we used a time step t = 3 s (Lesser et al., 2004). Each run simulated the morphological changes for 10 days under constant conditions during the simulation. To decrease the simulation time, we used a morphological scale factor equal to five. Aside from the Broad Creek model, the morphological factor value was chosen to avoid numerical instability.
The simulations were carried out by maintaining a constant wave height while everything else varied, for all simulated wave heights, then a constant wave period with everything else variable, and so on ( Table 2). We completed 144 simulations in total.

RESULTS
Our focus was on understanding the sediment exchange between the created marshes of living shorelines and adjacent SAV beds using the coupled Delft3D-SWAN model. In this section, we first compare the model results to field observations for validation, then present the hydrodynamic and morphodynamic results from the idealized model.

Model Validation
The estimation of the deposition rate made in the study by Bolton (2020) (performed along predefined transects) was carried out in the field from October 18, 2017 (HG and RU) to October 19, 2018 (OPP) for about 2 months. Then, annual deposition was extracted using lab techniques [see Bolton (2020) for further details]. Model validation was performed using the same timeframe as the sediment core samples. We ran two simulations for both years 2017 and 2018, from October 18 to November 3. The duration of the simulations was chosen in order to represent real conditions, but, at the same time, to have a reasonable computational time. The validation of the deposition rates was carried out along transects corresponding to 10 model cells, placed in correspondence with the measurement areas used in the study by Bolton (2020) (Figure 4). The obtained values were equal to the average obtained along the transects. The model predicted the deposition rates well in both years. Overall, deposition was underestimated for the three study sites but was still close to the real rate. Hatton Garden was more underestimated in 2017 but  was closer to the real value in 2018, while OPP and RU were close in both years ( Figure 5A). The SD showed significant variability of the deposition rate along the OPP transect and more contained variation along the HG and RU transects.
The underestimation of the deposition rate was likely attributable to the lack of waves in the modeling. Waves within the Chesapeake Bay are of the order of 20-30 cm under normal conditions and reach 1 m in height in extreme events (Lin et al., 2002), allowing and amplifying solid transport.
The water level extracted from the CBOFS was well-predicted by the model for both simulated years, showing high correlation coefficients (Figures 5B-E).

Hydrodynamic Results
Model results showed an average wave reduction between 10 and 80% in incoming wave heights (Figure 6). Wave damping was similar for the three different configurations. The addiction of the SAV did not induce particular changes in wave height reduction compared with the marsh only scenario, while the presence of rip-rap resulted in a slight but greater wave height reduction for higher wave heights ( Figure 6A). Wave damping was positively correlated with wave height following a linear law (Figure 6A), while shorter waves were dampened more than longer ones, which suffered a slight increase in height due to the shoaling effect ( Figure 6B). The slope of the seabed did not particularly affect wave reduction; in contrast, it was more important for shear stress. Higher slopes led to greater shear stress values (Supplementary Figure 3). The shear stress at the bottom was significantly reduced within the SAV bed, but then rapidly increases before being drastically dampened by the marsh on the shore ( Figure 7A). The shear stress at the beginning of the marsh (x = 200), for the configuration with only marsh on the shore, did not differ much from the scenario that also included the SAV, while the addition of riprap resulted in a reduction of shear stress values (Figures 7B,C) compared with both previous cases. Within the SAV bed, the marsh-only scenario produced the highest shear stress values (there was no SAV bed in the marsh-only scenario and, thus, less resistance to the flow), which were significantly reduced in the scenario that also included SAV. The addition of rip-rap resulted in shear stress values very close to the case with marsh and SAV (Figures 7D,E) and, thus, were lower than the marshonly scenario. Shear stress was positively correlated with slope, wave height, and period following a power law. Wave height, however, appeared to have a greater impact on shear stress than wave period, as evidenced by the higher values of the correlation coefficient R 2 (Figures 7D,E).

Morphodynamic Results
The sedimentation within the marsh was proportional to wave height and period for all simulated scenarios. The sedimentation for the marsh-only configuration was almost the same as the case with marsh and SAV, which did not significantly impact the deposition in the marsh. The marsh, SAV, and rip-rap scenario resulted in a slightly higher deposition for larger waves and slightly lower deposition for smaller waves (Figures 8A-C).
The deposition within the SAV bed varied according to the considered scenario. The case with the marsh only did not include any submerged vegetation and resulted in erosion within the location of the SAV bed, without any accumulation. The configuration with marsh and SAV, on the other hand, allowed the sedimentation inside the SAV bed to be proportional to wave height and inversely proportional to the period up to wave heights equal to 0.4 m. For wave heights equal to 0.5 m, deposition was also proportional to the period. Higher waves showed a decrease in deposition associated with greater shear stress values, which caused the erosion of the space in between marsh and SAV that we defined as "scarp, " also including part of the SAV bed ( Figure 8D). The scenario with marsh, SAV, and rip-rap showed less deposition in the SAV compared with the case with marsh and SAV, but it was greater than the marsh-only scenario. The accumulation within the SAV bed was proportional to wave height and period up to wave heights equal to 0.4 m. Higher wave heights and periods led to higher shear stress values and wave reflection due to rip-rap, which resulted in less deposition. However, wave heights equal to 1 m were able to transport the greatest amount of sediment when accompanied by longer periods, thanks to the sediment resuspension caused by the high shear stress values ( Figure 8E).
We calculated the erosion as the negative difference between the initial and final bed level (this included the marsh, the SAV, and the "scarp"). The erosion was proportional to both period and wave height for scenarios with marsh only and both marsh and SAV (Figures 8F,G). The scenario that also included the rip-rap reduced the erosion by an order of magnitude when compared with the other two cases. Erosion was proportional to wave height and period up to wave heights equal to 0.5 m. Higher wave heights (0.7 and 1 m) were capable of decreasing erosion since more deposition was associated with these conditions. Extreme events were capable of resuspension and transport the highest amount of sediments thanks to their high shear stress values, which allowed the replenishment of the shoreline with sediment, subsequently reducing erosion down to values close to zero ( Figure 8H).
On average, marsh deposition was similar for all simulated scenarios, wherein there was slightly higher deposition for the marshy-only scenario compared with the other two. The deposition in the SAV was reduced by the presence of rip-rap compared with the case with marsh and SAV. Erosion was slightly reduced in the marsh + SAV scenario compared with the marsh only, while it was drastically reduced in the case that also included rip-rap. The structure protected the marsh edge from being eroded, resulting in less erosion than the two cases without riprap but increased wave reflection, further resulting in the greater erosion of the SAV bed (Figure 8).
The increase in basin slope resulted in less deposition within the marsh, greater deposition in the SAV, and higher erosion compared with the gentler slope (Supplementary Table 1).
We then estimated the percentage of deposition in the marsh, in the SAV, and erosion for each of the three examined configurations with respect to the sum of deposition (marsh + SAV) and the absolute value of erosion of each scenario (Figure 9). We distinguished erosion according to whether it was a marsh, SAV or scarp (we defined scarp as the space in between the marsh edge and the SAV bed). The marsh-only scenario resulted in 31% deposition and 69% erosion ( Figure 9A). The 45% of the deposition in the marsh was caused by storms, while the erosion affected the marsh for 53% and the scarp for the remaining 47% ( Figure 9D). The marsh and SAV scenario reduced erosion and increased deposition compared with the marsh only case. It resulted in 28% deposition in the marsh, 15% deposition in the SAV, and 57% of the erosion (Figure 9B). The 44% of the deposition in the marsh was caused by storms, extreme events impacted the deposition in the SAV for only 34%, and erosion was divided into 60% marsh, 30% scarp, and 10% SAV ( Figure 9E). The marsh, SAV, and rip-rap scenario reduced erosion and increased deposition compared with the other two cases. It resulted in 50% deposition in the marsh, 16% deposition in the SAV, and 34% of the erosion (Figure 9C). The 60% of the deposition in the marsh was caused by storms, the extreme events impacted the deposition in the SAV by 27%, and erosion was divided into 70% scarp and 30% SAV (Figure 9F).
We calculated, starting from erosion and deposition results, the sediment balance that occurred between marsh and SAV among all simulated scenarios. The sediment balance was defined as deposition into marsh + deposition into SAV -erosion, and it was calculated using the control domain shown in Figure 10. The sediment budget quantification is shown in Figure 11 for the different wave heights. Results were averaged over the period. The sediment balance was negative for the configuration that only included the marsh on the shore, indicating a net loss of sediment that was higher for the increased slope (Figures 11A,D). The sediment balance with both vegetation communities drastically reduced sediment loss for low-energy and slope scenarios but still had a higher amount of sediment loss for higher energy and slope scenarios. The sediment budget was positive (close to zero) for almost all of the simulated cases with a low slope, except for the extreme wave heights, 0.7 and 1 m. The balance was negative at the higher slope even for waves equal to 0.3, 0.4, and 0.5 m (Figures 11B,E). The configuration with marsh, SAV, and rip-rap resulted in a higher sediment retainment for all wave conditions. The sediment budget was slightly positive under the action of ordinary external forcing. For higher wave heights, which caused the greatest loss of sediments in the other two configurations, the sediment balance was strongly positive for gentle slopes, while it resulted in less accentuation for greater slopes (Figures 11C,F). The effect of increasing the wave period was an increase in the loss of sediments (scenario with marsh only and both marsh and SAV) while allowing greater gains in the scenario that also included rip-rap. Extreme events were associated with larger wave heights, but also with larger periods (Supplementary Figure 4).
A meaningful summary of the effect of the three different configurations on sediment balance is shown in Table 3, where the budget values were averaged over the different wave heights and periods. The marsh + SAV configuration greatly reduced sediment loss compared with the marsh only scenario, and the aspect was made even more evident when rip-rap was also included.

Comparison With Previous Studies
Our numerical experiment, which couples Delft3D-SWAN, revealed the pivotal role of sediment interchange between the created saltmarsh of a living shoreline and its adjacent SAV. Simulation results demonstrated the efficacy of the vegetated shoreline in absorbing wave motion, thereby reducing incoming wave height by up to 80% according to a study by Manis et al. (2015), which found a wave energy reduction by Crassostrea virginica (eastern oyster) and Spartina alterniflora up to 67% in a wave tank experiment. Longer waves undergo an increase in height due to the shoaling effect, and according to a study by Battjes et al. (2004), they are nearly fully reflected at the shoreline, while higher-frequency components are subject to significant dissipation in narrow inshore zones including swash zones. Our model neglected wave reflection, so wave energy was all dissipated at the coast, which may result in greater erosion caused by longer waves. The model lacked the sensitivity to capture the SAV effect on wave height attenuation, which was mainly absorbed by the marsh, as well as by the friction at the bottom and wave breaking. However, previous studies have shown that submerged vegetation can attenuate wave height when vegetation height is comparable to water depth (Ward et al., 1984;Fonseca and Cahalan, 1992) and when wave orbital velocities and the seagrass canopy interact with each other (Chen et al., 2007). However, a study by Chen et al. (2007) highlighted that, to adequately assess the effect of SAV in coastal protection, its spatial and seasonal variability should be considered, together with the random variability of extreme events that mostly shape the coast and the spectral or directional distributions of wave energy. Our model kept the same wave direction orthogonal to the coastline for all simulations and the same vegetation density and height; thus, our results are limited to our study cases. Nevertheless, our model provided insight into the fundamental physical forces and principles underlying sediment transport and fate.
The higher shear stress values within the two vegetations were associated with longer periods, greater wave heights and basin slope, as found by Vona et al. (2020) regarding marshes and breakwaters interaction.
The modeling results showed the key role of waves in sedimentation and erosion. Extreme waves, with greater heights and periods, are capable of transporting the greatest amount of sediment within the marsh, but also have the greatest erosive power at the same time (Figure 8). Extreme events are also harmful to the SAV. The deposition was allowed more as the wave height increases, but storm conditions (wave height equal 0.7 and 1 m) erode the scarp mostly involving part of the SAV bed. However, the presence of rip-rap allowed deposition in the marsh as it greatly reduced erosion in case of extreme events, allowing the vegetation to obtain sediment supplies without suffering any losses. On the other hand, the artificial structure resulted in less deposition in the marsh in the case of normal external forcing. Studies in the literature have shown a similar behavior in waves when allowing sedimentation into the marsh (Castagno et al., 2018;Duvall et al., 2019;Nardin et al., 2020;Vona et al., 2020), also underlining the important role of other factors such as the alongshore current and the tide role in sediment retention within shallow coastal bays. On average, riprap resulted in a marked decrease in erosion. The marsh edge was protected by the structure and did not suffer any loss, but the rip-rap reflected incoming waves more, causing more erosion in the scarp and the SAV bed, which resulted in less deposition (Supplementary Table 1). Our model did not take into account the permeability of the break walls constituting rip-rap, since the structure was implemented in the model as an integral part of the bathymetry with a non-erodible bottom. The dissipative effects of the structure on transmitted waves were calculated using the wave transformation implemented in SWAN (wave-wave interactions, wave refraction, and wave dissipation by bottom friction and wave breaking), which did not take the effect of porosity into account. However, a study by Safak et al. (2020) investigated the effect of wave transmission through living shoreline break walls, finding that well-engineered semi-porous living shorelines act as buffers against human-mediated boat traffic and waves. They further underlined the benevolent aspects of coupling between natural and artificial solutions when studied correctly.
Our results highlighted the crucial role of SAV in coastal protection, as it reduces sediment export from the bay and allows sediment retention. This is in agreement with Donatelli et al. (2018), who highlighted the benevolent ability of seagrasses to the increase sediment storage capacity within shallow coastal bays. On the other hand, a decrease in seagrasses reduced the ability of the system to retain sediments, as our model shows when only the marsh occupies the shoreline, always determining a negative sedimentary balance. Donatelli et al.

Marsh
Marsh + SAV Marsh + SAV +Rip-rap Slope 0.8% −2.15 × 10 3 −8.9 × 10 2 9.9 × 10 2 Slope 1% −3.4 × 10 3 −1.8 × 10 3 3.7 × 10 2 (2018) also pointed out that the presence of seagrasses reduces the suspended sediment concentration available for the marsh, with this aspect being slightly captured by our model as the deposition in the marsh was slightly higher in the absence of the SAV bed (Supplementary Table 1). However, a study by Chen et al. (2007) pointed out that, in order to fully understand the sediment retention mechanisms by the SAV, different sediment concentration inputs and transport within the bed due to currents should be taken into account. Our model maintained a constant inlet concentration of 0.1 kg/m 3 . A higher concentration would lead to more sedimentation in the marsh, as evidenced in the study by Vona et al. (2020), while deposition in the SAV would require further studies. The presence of a non-erodible structure such as rip-rap slightly affected the sediment budget in the case of normal external forcing while helping stabilize the shoreline under extreme events. In storm conditions, solids transport in coastal wetlands was magnified, which implies the presence of riprap can be crucial in avoiding greater damage while allowing extreme events to replenish marshes and vegetated shorelines with sediments, as found in the studies by Castagno et al. (2018) and Vona et al. (2020).
The modeling results showed that hybrid infrastructures harnessed the benefits of both natural and built solutions to improve shoreline resilience. Moreover, by including vegetated features, hybrid systems are not likely to lose their effectiveness over time due to SLR, as in the case of submerged breakwaters for instance. This is because they can adapt themselves to SLR as long as they receive the right sediment replenishment, as shown in the study by Sutton-Grier et al. (2015).
Another important aspect to consider to improve the living shoreline configuration and make it more natural, while maintaining its efficiency, is the possibility of replacing the break walls constituting the rip-raps with oyster reefs. The advantage of integrating in-water infrastructure with oysters is that these reefs have the ability to grow with SLR (Rodriguez et al., 2014;Ridge et al., 2017), providing greater guarantees in terms of coastal protection even after the gray structure may become ineffective.
Interactions between vegetation species will play a key role in determining coastal wetlands responding to SLR and extreme events, and it will be crucial to keep understanding and monitoring feedback between different plant communities to better safeguard our coastal environments.
However, the evolution of vertical and horizontal marshes over time is a topic that still requires future research. Vertically, marshes need to accumulate sediments to contrast SLR, while horizontally, they need landward expansion to compensate for erosion (Fagherazzi et al., 2020). The dynamics of marsh evolution will require future studies to provide us with tools to safeguard our coasts.

Model Limitations
Delft3D-SWAN is a three dimensional, depth-averaged software for hydrodynamic computation. Hydraulic roughness due to vegetation is modeled, for rigid vegetation, by Baptist equations (Baptist, 2005), while for flexible vegetation, Delft3D assumes a greater degree of roughness (Lera et al., 2019).
However, models did provide useful insights on sediment transport and fate. In fact, we were able to model the mass balance of sediment distribution by coupling the hydrodynamic module with the vegetation model in the study of Baptist, with the possibility of adding different sediment characteristics Vona et al., 2020).
Baptist's equation has been widely tested with field and laboratory experiments, such as in the Allier River (France), the Volga River (Russia), and the Rhine River (Netherlands), with natural and artificial vegetation, coming to the conclusion that the model predicted sedimentation differences caused by vegetation well. Many experiments, compared with the results in the study of Baptist, gave back comparable outcomes. Recently, the study of Crosato and Saleh (2011) provided another validation of the Baptist equation with field observations applied to the Allier River in France. Moreover, the study of Baptist et al. (2007) used the results of the depth-averaged kappaepsilon turbulence model, which accounts for vegetation in a genetic programming framework, to obtain an expression for roughness in the presence of vegetation, using a variety of input parameters to find the dimensionally consistent, symbolic equation.

Guidelines for "Living Shoreline" Implementation and Success
Coastal communities face ever harder challenges due to the rise of extreme events and SLR. Coastal stabilization solutions do not only require the implementation of engineering techniques, such as breakwaters or bulkheads, since new stabilization options, such as the living shoreline, can reduce erosion and man-made structures while providing ecosystem services such as food production, nutrient removal, and water quality improvement (NOAA, 2015).
However, the success of vegetation planting for coastal protection is not always guaranteed, as it depends on various factors that can affect plant resilience. One of the most important factors is wave climate. Wave climate is influenced by factors such as fetch, wind speed and duration, water depth, basin slope, and the orientation of the site. When a shoreline is facing the storm wind direction, it is more likely to be damaged, while a shallow water environment with a gentle slope is more capable of reducing wave energy and protecting plant species. Other important factors for the living shoreline survival are constant sediment supply, given by periodic tidal flooding, nutrient supply and salinity, and many others that can be found in the study of Broome et al. (1992). Given the dependence of the planting success on the geographical site, the study results highlighted how the ideal conditions are found for living shoreline installation, as the simultaneous presence of marsh and SAV greatly helps to reduce the sediment loss outgoing coastal wetlands. This improves the resilience and ensures the better stabilization of the shore, while riprap considerably helps with storm conditions (Figure 12). Our work can give good indications on the implementation of coastal recovery and maintenance plans in shallow coastal bays, given wave climate and bathymetric conditions, through natural solutions that look at both the protection and the ecology of coastal environments.

CONCLUSION
Understanding sediment transport dynamics is crucial for coastal protection. Our study provided references for decision makers working in coastal wetland restoration. In this study, we investigated the sediment interchange between saltmarshes and SAV to better quantify the synergy between the two plant communities. The SAV certainly helped enhance shoreline resilience, reducing sediment loss by retaining the sediments outgoing the marsh and allowing deposition. The presence of non-erodible structures such as rip-rap, which is often an integral part of the living shoreline, protects against extreme events and allows sediment replenishment inside the marsh, especially under storm conditions. The delicate balance between the supply and loss of sediments is crucial in coastal wetland restoration. Artificial structures make it possible to avoid greater damage in extreme events, which also has a significant economic impact. Living shorelines certainly offer a valid and green solution to coastal protection, which, if coupled with man-made or natural structures, such as oyster reefs rather than rip-rap, can strengthen the resilience and the vitality of the coast.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.