The modeling of marsh erosion due to lateral wave thrust is implemented in the open-source Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) numerical modeling system (Warner et al., 2010). The COAWST modeling framework couples the ROMS (Regional Ocean Modeling System) model for hydrodynamics with a wave model - SWAN (Simulating Waves Nearshore) via the Model Coupling Toolkit (MCT) (Warner et al., 2008). ROMS is a three-dimensional, free-surface, finite-difference, terrain-following model that solves the Reynolds-Averaged Navier-Stokes equations using the hydrostatic and Boussinesq assumptions (Haidvogel et al., 2008). SWAN (Simulating Waves Nearshore) is a third-generation spectral wave model based on the action balance equation (Booij et al., 1999). After a user-defined number of time steps, there is an exchange of water level and depth-averaged velocities from ROMS to SWAN and wave fields from SWAN to ROMS. The Community Sediment Transport Modeling System (CSTMS) model accounts for the 3-D transport of sediment using the bedload and suspended-sediment components. The bedload mass can be updated using a variety of parameterizations that require bed shear stress based on current and wave forcing from the bottom cell. The suspended-sediment is transported by solving an advection-diffusion equation which accounts for a source/sink term that leads to a vertical exchange or settling with the bed. The details of these methods are explained in Warner et al. (2008). The model can represent any number of user defined sediment classes divided into cohesive and non-cohesive types. The amount of sediment stored in the bed is determined through the user-defined properties of each sediment class and sediment bed layer. These properties include bed thickness, sediment density, bed thickness, and bed porosity. In addition to the sediment transport model, the COAWST modeling framework can also account for the change in current and wave dynamics due to the presence of vegetation (Beudin et al., 2017).

Beudin et al. (2017) implemented the physical effects of vegetation in a vertically varying water column through momentum extraction, vertical mixing, and wave dissipation. This allows for the modeling of the impact of marsh vegetation in preventing marsh surface erosion through wave energy dampening and sediment trapping (Möller, 2006; Le Hir et al., 2007; Mariotti and Fagherazzi, 2010). In the current work, we implemented the method to account for marsh edge erosion through the action of lateral wave thrust (LWT) and the method to account for vertical growth of marsh complex due to biomass production. Figure 1B shows the flowchart of the coupled modeling framework of COAWST with the addition of routines for modeling marsh dynamics.

The presence of marsh is defined through a user-defined mask provided to the model as input. A value of 1 corresponds to marsh cells, while a value of 0 is associated with non-marsh regions. Note that this masking operation is independent of the wetting and drying masking framework of COAWST (Warner et al., 2013). The marsh masking can change from 1 to 0 at a given cell once a given amount of user defined marsh retreat due to lateral erosion occurs in the model; this implies that the cell is no longer a part of the marsh. Once the marsh cells are converted to non-marsh cells, they are retained as non-marsh cells.

Once the user specifies the initial marsh mask, the boundary of marsh/non-marsh regions is identified during model runtime. The wave thrust per unit width is calculated by taking the vertical integral of the dynamic wave pressure. The wave thrust is divided into above and below mean sea level components based on the formulations described in the Department of the Army, Waterways Experiment Station, Corps of Engineers (1984). This formulation to compute wave thrust on marsh has been widely utilized in earlier works (Tonelli et al., 2010; Francalanci et al., 2013; Bendoni et al., 2014; Leonardi and Fagherazzi, 2014). Next, we describe the formulation of lateral wave thrust in the model.

First, the above mean sea level component that accounts for the hydrostatic pressure from wind waves (LWT_ASL) is computed and is defined as:

where ρis the density of water, g is the acceleration due to gravity, and H_s is the significant wave height. Second, the below sea level component, LWT_BSL, accounts for the dynamic pressure of wind waves by including the effect of changing water level and is defined as:

where K_p is the pressure-response factor due to water particle acceleration under the effect of wind waves and is calculated as:

whereh is the elevation of the marsh platform at the edge of the scarp with respect to mean sea level, ζ is the water level, and k == 2⁢πLw⁢a⁢v⁢e, where L_wave is the wave length. Next, the thrust from the two components can be added (Eqs 1, 2) to give the total thrust due to wave attack.

Based on the direction of the waves and grid orientation, the fraction of total thrust that is normal to the marsh cell faced is determined and used in subsequent calculations. All nearest neighbors of a marsh cell can provide LWT (with a maximum of four neighboring cells in the structured grid approach) and then the total wave thrust from all neighboring cells is added to give a total wave thrust on the marsh cell (Supplementary Figure 1 shows the application on LWT on marsh cells).

Next, the effect of water overtopping the marsh on wave thrust is included. Tonelli et al. (2010) showed that the wave thrust reached a maximum when the water level was co-located with the marsh scarp elevation and reduced when waves overtopped the marsh scarp. To account for the change in wave thrust based on water level, various studies have formed different parameterizations (Tonelli et al., 2010; Leonardi et al., 2016b). In the current model implementation, the wave thrust is reduced exponentially as water level increases above the marsh scarp. After modifying the wave thrust based on water level, the resulting wave thrust magnitudes from all the cell faces, i.e., for marsh cells that have multiple edges exposed to the estuary, are summed to obtain a total thrust at each cell center.

The mass of sediment released from marsh cells (m_LWT–exportin kg) depends on the lateral wave thrust (LWT in kN/m), marsh erodibility coefficient (k_marshin s/m), grid size (dx in m), and time step size (dt in s).

Based on these factors, the erosion of sediment occurs through a change in bed mass by taking the sediment out from the marsh cell and adding it to the adjacent cell’s bed mass. In the case of a marsh cell surrounded by multiple neighboring grid cells that can provide wave thrust, the sediment change occurs in accordance to Eq. 5. The changes in bed mass modify the bathymetry of the domain leading to a change in bed morphology. Note that the sediment stored in the marsh is assumed to be on the top bed layer, and the model does not account for erosion of material from deeper seabed layers.

As mentioned in the previous section, the action of LWT results in a loss of sediment from the marsh. Once enough sediment is exported from the marsh to account for a user-defined reduction in the scarp height, the marsh cell is considered to undergo a horizontal retreat that is equal to the width of the marsh cell. The marsh cell then converts into an open water cell. Specifically, the marsh mask is changed from 1 to 0 for this cell. Note the reduction in scarp height is based on the absolute change of bed thickness. As the marsh cells convert to open water cells, the vegetation in the marsh cells also suffers a dieback and this is simulated in the model by setting the vegetation biomass to zero. Vegetation regrowth is not allowed in marsh cells converted to open water cells. The sediment in the open water cell that was previously a marsh cell can get transported under the action of hydrodynamics.

The broadly used (Mariotti and Carr, 2014; Kirwan et al., 2016; Carr et al., 2018, 2020; Mariotti, 2020) growth model formulation is adapted from Morris et al. (2002) following the formulation of Kirwan and Murray (2007) where vegetation rapidly adjusts to changes in elevation and as such biomass can be calculated as a function of marsh cell depth below a Mean High Water (MHW). In these two studies, the biomass productivity formulations were based on measurements from S. alterniflora. Biomass productivity (B_peak)is based on a parabolic biomass curve where the upper (D_max)and lower limits (D_min) of the parabola are a function of MHW and is defined as:

where B_max is the optimal biomass in kg m^–2yr^–1 that is a user input where

and where MTR is the mean tidal range and is assumed to be MTR = 2MHW. cff in the denominator of Eq. 6 is a scaling factor that does not allow the value of B_peak to exceed a maximum value of B_max set as a user input of 2.5 kg m^–2yr^–1 and is defined as:

Meanhighwater is calculated internally as a moving average over user defined days by keeping a track of the maximum water level in a day. The upper and lower limits correspond to reference depths where the macrophyte survives and leads to accretion of biomass (McKee and Patrick, 1988). The integrated per year amount of below ground biomass (AMC) corresponding to 180 days of growth in kg m^–2yr^–1 is calculated as:

where nu_gp is the fraction of below ground biomass set to be an input of 0.0138 day^–1 and tdays_growth is set to be 180 days. tdays_growth represents a growing period for marsh vegetation within the marshes in Barnegat Bay. The effective biomassR_ref in kg m^–2yr^–1 after accounting for recalcitrant carbon is calculated as:

where chi_ref is the fraction of recalcitrant Carbon set to be an input of 0.158. The model choice of nu_gp, tdays_growth, and chi_ref are all based on Mudd et al. (2009). The rate of vertical growth due to biomass production over marsh cells (m/year) is calculated from R_ref using:

where m_bulk–dens is the bulk density of marsh organic matter. The vertical growth rate is used to calculate vertical accretion and is then converted to a change in bed mass in marsh cells.

The numerical experiments of Tonelli et al. (2010) calculate maximum wave thrust for three different wave heights and provides a dataset to verify the implementation within the current model. The setup of the domain involves a basin of 10 by 1 km rectangular domain leading to a grid resolution of 100 by 50 m with a depth of 1.3 m. The plan view of the domain is shown in Figure 2A and the variation in bathymetry along the channel is shown in Figure 2B. The marsh complex is located at an elevation of 0.1 m above mean sea level. The elevation of 0.1 above MSL is based on the numerical experiments of Tonelli et al. (2010) where the mean water level of the numerical experiments is set to 1.4 m above the bottom of the domain. The model is forced by oscillating the water level on the northern edge with a tidal amplitude of 0.5 m and a period of 12 h. The waves are imposed in the entire domain to have a constant wave height and wavelength. The waves impinge at a 0 degrees angle with respect to the northern boundary with a period of 2 s. Because the waves are constant in the domain, this simulation remains uncoupled (no wave model). The constant waves in the domain are generated by providing bulk wave statistics in an input file to match the wave conditions of Tonelli et al. (2010). The ROMS barotropic and baroclinic time steps are 20 and 1 s, respectively. The wave height is varied to find the variation of maximum wave thrust with water level and compared with the results from Tonelli et al. (2010). The model setup is run for two tidal cycles allowing computation of the influence of varying water levels on the maximum lateral wave thrust. Three different scenarios of wave height are simulated by forcing the northern boundary with heights of 0.2, 0.3, and 0.4 m.

An idealized test case with higher spatial resolution is developed to test the model’s ability to simulate marsh lateral retreat on the order of 1 m. Figure 2C shows the plan view of the model domain that is 600 m long and 650 m wide with a grid resolution of 1 and 25 m in the cross-shore and along-shore directions, respectively. A grid resolution of 1 m allows for modeling a realistic case of retreat of marsh cells that occurs from a monthly to annual time scale while also simulating short timescale dynamics with a timestep of 1 s. The initial bathymetry consists of a maximum depth of 1.6 m corresponding to the seaward side (eastern boundary) that tapers to 0.2 m depth spanning over 353 m of the domain. Beyond that, a fixed elevation of 0.35 m describes the start of the marsh complex. The variation in bathymetry along the length of the channel at a cross-section (y = 325 m) is shown in Figure 2D. The marsh complex in the setup is vegetated and one sediment class is used for the entire domain with marsh vegetation and sediment properties presented in Supplementary Table 2. The model is forced by oscillating the water level on the eastern edge with a tidal amplitude of 0.5 m and a period of 12 h. Waves are also imposed on the eastern edge with a height of 0.3 m, directed to the west (90 degrees) with a period of 6.28 s. The northern and southern boundaries of the domain are closed. The bottom boundary layer roughness is increased by the presence of waves that produce enhanced drag on the mean flow (Ganju and Sherwood, 2010). The turbulence model selected is the k−ε scheme (Warner et al., 2005). The complete list of model parameters is shown in Supplementary Table 2.

To initialize the model, the simulation is setup without the action of LWT on the marsh domain. This allows for erosion to occur from the elevated platform due to only wave-current stresses (WCS). The model is run until the marsh platform topography has stabilized. The simulations are then restarted with the introduction of marsh mask and inclusion of LWT routines. After that the marsh edge erodes under the combined effect of WCS and LWT. These simulations are continued until the marsh edge cell retreats and a second marsh edge cell is formed using the methods described above.

We test the combined effects of dynamic sediment transport, lateral wave thrust induced erosion, and vertical growth due to biomass productivity using two different simulations involving marsh complexes within Forsythe National Wildlife Refuge and the Barnegat Bay (BB) estuary. The COAWST BB model was developed and applied by Defne et al. (2019). We use the same setup of the BB domain (Figure 3A) and simulated the model grid for the period of May, 2015. The tidal forcing is generated from the ADCIRC product (Szpilka et al., 2016) and wind forcing is generated from the 3-hourly North American Mesoscale Forecast System product (NCEP NAM, 2020). Three estuarine sediment types and one marsh sediment class were utilized following the study of Defne et al. (2019). The marsh and submerged aquatic vegetation (SAV) presence along with their properties in the estuary was determined using the model development of Donatelli et al. (2018).

We used the model output from the BB simulation to provide the hydrodynamic (Supplementary Figure 4) and sediment forcing (water level, barotropic and baroclinic velocities, waves, and suspended sediment) to provide hourly boundary forcing conditions for the two creek grids. The grid dimensions for the two marsh complexes are 132 × 388 points and 599 × 539 points, at 5-m resolution with 7 vertical sigma layers, for Reedy Creek and Dinner Creek, respectively (Figures 3B,C). The northern, southern and eastern boundaries of the two domains remain open to BB domain forcing, while the western boundaries are closed. Initial model bathymetry is obtained from the CoNED topo bathymetric model for New Jersey and Delaware (OCM Partners, 2021). The bathymetry in the channel is modified to be 1 m in the Reedy creek domain and 2.3 m in the Dinner creek domain to match the observed bathymetry of Suttles et al. (2016). The marsh elevation from CoNED for the RC domain is adjusted by 0.1 m to field measurements (real time kinematic global positioning, Quirk, 2016) and comparison to local mean sea level. A similar elevation shift (∼0.1 m) is seen in the local mean sea level computed from a decade long water level time series at the nearby Mantoloking station. For the DC domain, we retained the marsh elevation as provided by the CoNED model. For each domain, the cells with marsh mask are based on the updated National Wetland Inventory delineation for BB (Defne and Ganju, 2016) and defined in the model as all the cells covering the vegetated areas or cells under a depth of 0.1 m. The same area is used to specify marsh vegetation properties (Supplementary Table 3). The suspended marsh sediment class in BB domain is exchanged with the creek domains as part of the open boundary forcing. In addition to this, two more marsh sediment classes, organic and inorganic, were introduced in the creek simulations. In the simulations, the organic marsh class can accumulate only through vertical growth processes. This allows for distinguishing the sources of sediment on the marsh complex that could be derived from either organic (first) or inorganic (second) class of sediment or the inorganic (third) class available from the BB domain boundary. Supplementary Table 3 specifies the properties of the sediment classes, marsh/estuarine bed, and marsh vegetation used during the creek simulations. At the boundary of RC domain, the maximum and mean suspended sediment concentration (SSC) corresponding to the finest estuarine sediment class are 15 and 3 mg/L, respectively, while the maximum and mean SSC for the coarsest estuarine sediment are 0.2 and 0.002 mg/L. At the boundary of DC domain, the maximum and mean suspended sediment concentration (SSC) corresponding to the finest estuarine sediment class are 38.6 and 13 mg/L, respectively, while the maximum and mean SSC for the coarsest estuarine sediment are 3.9 and 0.45 mg/L. In the creek simulations, the ROMS time step is 1 s, while the SWAN time step and the coupling interval between ROMS and SWAN are 20 min. The friction exerted on the flow by the bed is calculated using the bottom boundary layer formulation (Warner et al., 2008). The bottom boundary layer roughness is increased by the presence of waves that enhance drag on the mean flow (Ganju and Sherwood, 2010). The vegetative drag coefficients (c_D) in the flow model and the wave model are set to 1. The bed roughness is set to z₀ = 0.015 m, which corresponds to a mixture of silt and sand (Soulsby, 1997). The turbulence model selected is the κ − ε scheme (Warner et al., 2005).

Under the influence of varying water levels in the model setup described in Section “Idealized Domain to Formulate and Verify Lateral Wave Thrust,” the maximum LWT is computed in three different scenarios corresponding to wave heights of 0.2, 0.3, and 0.4 m (Figure 4). The use of constant wave forcing in the current approach ensures that we can simulate wave conditions similar to the numerical experiments of Tonelli et al. (2010). LWT increases with increasing wave heights. As the water level goes above the marsh scarp, the LWT starts to reduce exponentially because water runs on top of the marsh scarp. LWT peaks when the water level elevation is equal to the scarp height. As the water level drops, the integrated LWT starts to decrease over the marsh scarp. The pattern of LWT variation with depth above bottom is consistent with the numerical experiments from Tonelli et al. (2010). However, because the present approach relies on bulk wave statistics to parameterize LWT (section “Computing LWT”) the degree of non-linearity is much lower in LWT variation. On the other hand, the numerical experiments of Tonelli et al. (2010) capture the non-linear wave dynamics in their modeling methodology and calculate LWT directly by integrating pressure forces over the marsh scarp.

Initially, model setup is simulated without the action of LWT while an elevated platform starts to erode under the action of combined WCS until the bathymetry at the edge of the elevated platform is stable. The transect at x = 248 m shown in Figure 5A corresponds to the cell that dropped the most in elevation only due to WCS (marsh edge pre-LWT) and provides the topography at the end of 351 days, where the marsh edge pre-LWT had a change in elevation from 0.35 to 0.02 m. The cells at or below the depth of 0.02 m are considered open water cells and the remaining cells in the cross-channel direction (east to west direction) that still maintain an elevation greater 0.02 m in the domain are masked to be marsh cells. At this time, the edge of the marsh complex is at an elevation of 0.27 m starting from the marsh edge (x = 247 m, referred to as “pre-LWT” in Figure 5A). The user defined scarp height to erode the edge due to the lateral wave action is also set to this elevation (i.e., 0.27 m with the marsh edge located at x = 247 m). After setting the marsh mask for the marsh complex, the simulation is restarted with the effect of LWT included.

Following the WCS-only run, the model is resumed with the LWT formulations. For the same transect described in the previous paragraph (marsh edge 1 at x = 247 m), the marsh edge drops to an elevation of 0.13 m after 216 days. After 216-days of model run, the marsh edge (x = 247 m) drops to an elevation of 0.13 m. The combined processes of WCS and LWT change the bed thickness to the user defined scarp height of 0.27 m. Based on the retreat formulation (section “Evolution of Marsh Mask and Computing Lateral Marsh Retreat”), the marsh edge cell (marsh edge 1 post-LWT at x = 247 m) converts to an open water cell, forming a unvegetated mudflat adjacent to the marsh (Figure 5B). Due to the retreat of the marsh edge cell at x = 247 m, a new marsh edge cell is created at x = 246 m (referred to as “marsh edge 2 post-LWT” at x = 246 m).

A change in the bed thickness of the sediment bed in the cross-channel direction during the erosion of marsh edge 1 post-LWT is shown in Figure 6A. The lowering of marsh complex elevation results in a lower vertically integrated wave thrust and therefore the total wave thrust decreases over time (Figure 6B). After marsh edge 1 post-LWT gets converted into an open water cell, the new marsh edge (marsh edge 2 post-LWT) is influenced by LWT action. A similar pattern is observed in the change in daily flux out of the marsh edge cells (Figure 6C), where the decrease in wave thrust with lowering elevation results in lower flux out of the marsh edge cells.

We tested the combined dynamics of lateral wave thrust (LWT)-based erosion and vertical growth due to biomass productivity using two different simulations involving creeks in the Barnegat Bay estuary. Figures 7A,B show the spatial pattern of LWT that is summed over the 31-day long simulation period for each of the two domains. In both domains, the cells closer to the open bay experience a higher LWT compared to the cells inland. The spatial variation in LWT is explored by looking at two points in the RC domain (Figure 8A). The first point close to the open bay is referred to as “Bay point” and the second point inland of the RC channel is referred to as “Inland point.” Figure 8B shows the significant wave height (H_s) generated by winds at the points referenced in Figure 8A. The root-mean-square (rms) H_s over the 31-day period for the bay point is 0.03 m while rms H_s for inland point is nearly zero. Other than H_s, the wave direction impacts the formulation of LWT as it only allows for the normal component of waves (after accounting the orientation of the grid) to provide wave thrust from adjacent wet cells to the dry marsh edge cells. The bay point receives a higher integrated LWT over the 31-day period with a rms LWT of 0.034 kN/m while rms LWT for inland point is nearly zero, highlighted in Figure 8C.

Next, we analyze the vertical growth rate (VGR) of organic sediment due to biomass productivity at the end of the simulation period for the two creek domains in Figure 9. The difference in VGR between the two domains is relatively small. The VGR increases north of the DC channel more than in the southern end while remaining high along the channel for the RC domain.

We use two points (inland and at the bay end) from the RC domain to understand the reasons for spatially varying VGR (Figure 10A) and the reference datums that describe the upper (D_min) and lower (D_max) limits of the parabolic growth curve along with the elevation of the two chosen points (Figure 10B). The rolling average to get D_min and D_max employs a user-defined period of 7 days. The bay point is at a lower elevation around 2 cm while the inland marsh point is 12 cm above datum. Following the parabolic formulation, the highest amount of growth occurs when the cell elevation is midway between the upper and lower limits of the parabola and decreases away from the middle datum. In Figure 10B, the elevation of the bay cell is close to D_max while the elevation of inland cell is located close to the midway point between D_min and D_max, providing an optimum elevation for biomass productivity within the salt marshes. Consequently, at the end of the simulation period, the inland cell grows at a 3.6 mm/year rate, higher than the 2.1 mm/year growth rate for the bay point. It should be noted that the VGR per year is computed based on a monthly calculation and would vary once the effects of seasonality on biomass productivity and on hydroperiod parameters are considered in a year long simulation. To compare the LWT and VGR pattern between the two creeks, we extracted the points along the channel from the two domains and categorized them into eastern and western points for the two domains (Figures 11A,B). For both the domains, the orientation of the waves results in a higher integrated LWT at the eastern end of the channels (Figures 11C,D). Due to the open mouth of the DC channel, the maximum integrated LWT over the 31-day period is higher for the eastern end of the DC domain (74.1 N/m) compared to the eastern end of RC domain (60.4 N/m). Between the western ends, integrated LWT is higher for the RC domain (16.4 N/m) compared to the DC marsh complex (12.2 N/m). In addition, the effect of LWT extended to a distance of 1000 m for the DC domain compared to 200 m for the RC domain. A comparison of VGR for the two domains made by analyzing the VGR variation with elevation show the mean elevation for the eastern and western ends of DC domain is higher (12.6 and 14.7 cm, respectively) compared to RC domain (11.9 and 13.9 cm, respectively). In addition, the tidal range for the DC domain is 70 cm while the RC domain is 20 cm. The combination of higher elevation marsh cells along with increased tidal range provides a larger envelope for the parabolic formulation growth curve for the DC domain compared to the RC domain (Figures 11E,F).

Next, we analyze the amount of sediment deposition from various processes in the two domains (Figure 12). This includes the amount of sediment deposition due to biomass productivity (Figures 12A,B), inorganic marsh sediment (Figures 12C,D) and inorganic estuarine sediment (Figures 12E,F) deposited from the process of tides and waves (hydrodynamic forcing) overtopping the marsh cells. In the RC domain the organic sediment is accumulated (Figure 12A) following the pattern of spatial changes in VGR due to biomass productivity observed in Figure 9A. The inorganic marsh sediment deposited from hydrodynamic forcing is limited to only a few cells (Figure 12C). The inorganic estuarine sediment deposited over the marsh is significantly higher in the southern end of the channel and decreases north of the RC domain (Figure 12E). Similar to the RC domain, the organic sediment accumulation follows the pattern of VGR for the DC domain (Figure 12B). It is observed that a larger number of cells accrete inorganic marsh sediment due to the hydrodynamic forcing (Figure 12D) and there is a larger deposition of inorganic estuarine sediment toward the northern end of the channel in the case of DC domain (Figure 12F).

Numerical models used to assess salt marsh trajectory typically incorporate a subset of the relevant physical processes (Fagherazzi et al., 2012; Ganju, 2019). More robust models must be developed to simultaneously account for lateral and vertical processes, in addition to dynamic sediment transport within the surrounding estuary, tidal channels, and intertidal flats. For example, vertical biogenic models (Morris et al., 2002; Alizad et al., 2016) alone will neglect the import of estuarine sediment along tidal channels that aid in landward deposition, while abiotic sediment transport models (Temmerman et al., 2005; Donatelli et al., 2018; Defne et al., 2019) neglect vertical response of the marsh to tidal levels. The modeling framework developed and applied here considers the dominant dynamic, lateral, and vertical processes in a coupled geomorphic context. This is necessary to properly evaluate marsh response to future climatic (morphodynamic) changes.

Within this modeling framework, there are three distinct sources of sediment that contribute to marsh development and deterioration. In addition to the process of sediment transport induced from wave-current stress, the present framework incorporates (1) the processes to account for the organic material accumulated due to biomass production and (2) the lateral wave thrust action leading to erosion of sediment from the marsh onto the estuary (Figure 1A). Table 1 shows the total mass of sediment deposited over all the marsh grid cells in the Reedy and Dinner Creek simulations. The largest source of sediment to the marsh in terms of mass in both creek simulations is the inorganic estuarine sediment (61% of deposition over the RC marsh and 94% of deposition over the DC marsh). Meanwhile, the total mass of inorganic marsh sediment (from marsh edge erosion) represents 25% of the total sediment accumulation over the marsh in RC and 3.7% in DC. The total mass accumulated over the Dinner Creek marsh is 45 times larger than the mass accumulated over Reedy Creek. This is consistent with the more energetic tidal currents and wave climate present in the Dinner Creek domain.

Figure 13 shows the extent of coverage of marsh grid cells from the three sources of sediment in the two creek domains More than 90% of the area of deposition in the marsh complex in both domains is dominated by the organic marsh sediment (Figure 13). In the entire RC marsh complex, the coverage of inorganic marsh sediment is over 0.4% of the marsh while the area covered in the DC marsh complex is 0.97%. In the DC domain, the inorganic estuarine sediment covers (5.7%) of the area compared to 5.3% in the RC domain (Figure 13). The comparison between the two creek domains shows that a larger area of marsh complex in the DC domain gets affected by the inorganic marsh and estuarine sediment compared to the RC domain. This can be attributed to the higher LWT and higher tidal forcing in the DC domain that leads to a higher availability of marsh and estuarine-derived inorganic sediment transport. The larger tidal range, and higher velocities therefore in the DC domain is also responsible for an increased accretion in the north end of the DC channel, where suspended-sediment concentrations are large due to advection (Figure 13B). When only the marsh areas dominated by the two inorganic sources of sediment (marsh and estuarine) are considered, the estuarine sediment dominates over 75% of those locations in RC domain. In contrast, 50.4% of the marsh cells on the DC marsh complex are dominated by inorganic sediment have a predominant marsh-edge erosion sediment source. These results mirror the higher tidal and wave energy in the southern part of Barnegat Bay, where sediment fluxes are larger leading to higher sediment load (Ganju et al., 2014) and likely increasing estuarine sediment transport onto the marsh platform. Organic accumulation dominates over most of the interior marsh plains, due to limited ability of sediment to be advected there given low or negligible water velocities and settling of sediment (Coleman et al., 2020). The large area of accretion due to biomass production shows that this process is the principal mechanism through which marsh complexes vertically accrete similar to the observations (Neubauer, 2008; Morris et al., 2016). The inorganic sediment (both estuarine and marsh) gets redeposited close to the marsh edge with a reduced deposition in the interior of marsh replicating the aspects of previous field observations (Reed et al., 1999; Temmerman et al., 2003; Roner et al., 2016; Lacy et al., 2020; Smith et al., 2021). Nonetheless, these results suggest that deposition on eroded edges help maintain elevation as the marsh contracts at the seaward edge (Hopkinson et al., 2018). Neglecting this aspect will give unrealistic spatial patterns in evolution that are not consistent with conceptual geomorphic models (FitzGerald et al., 2018). This highlights the importance of accounting for dynamic sediment transport in the study of marsh morphological modeling. Similarly, the advantage of having multiple sediment classes along with the inclusion of modeling marsh accretion and LWT based erosion provides the capability to trace the source of sediment (organic or inorganic marsh or estuarine sediment) that is available for export through estuarine dynamics.

The current methodology has several limitations that are listed below:

Beyond the implementation presented in this article, ongoing work has included two additional functions that are available for a variety of applications. First, the growth of marsh vegetation based on peak biomass using relationships mentioned in D’Alpaos et al. (2006) to calculate the change in stem density, stem length, and diameter based on peak biomass, implementing the feedbacks between growth of biomass and marsh vegetation properties. A second functionality involves the conversion of non-marsh cells (open water cells) into marsh cells, thereby allowing for modeling of marsh lateral expansion over tidal flat areas and allowing marsh vegetation cover to increase in an area. The modification of marsh vegetation properties based on biomass and creation of marsh cells could be used to explore the role of vegetation properties in altering marsh dynamics or studying the relationships between sediment supply and marsh coverage. Using these two functionalities along with morphological acceleration, an ongoing model application involves the study of long-term marsh formation under varying conditions of RSLR and sediment supply.