Influence of Belowground Biomass on Barrier-Island Evolution During Storms: A Computational Parameter Study

To assist in isolating the influence of dune topography and belowground biomass, we based our simulations on idealized alongshore-uniform topography. Selection of idealized topographic cross-shore profiles were informed by characteristics of the barrier islands along Moriches Bay, New York, USA. Ten cross-shore profiles were created from topographic lidar data (USACE et al. 2022) and bathymetric ETOPO data (NOAA 2022) in the Moriches Bay region (Cheng 2015). The lidar topography displayed large variability alongshore. For example, dune elevations range from approximately 4 to 8 m above mean sea level (MSL), distances from the shoreline to the dune crest range from 160 to 260 m, and total subaerial barrier island widths range from 100 to 600 m. In selecting idealized topographic characteristics, we gave preference to conditions likely to overwash and inundate during storms.

Fig. 2 provides an example of the idealized initial barrier island profile. The submerged portions of the initial profile in all simulations was fixed and specified by power laws with an exponent of 0.67 (Dean equilibrium profile) and constants of 0.15 and 0.02, respectively, in the ocean and in the bay. The initial subaerial profile in all simulations was characterized by a fixed cross-shore width of $x_b=500$ m and a fixed subaerial sediment volume of $V=688$ m³/m. The initial subaerial beach and back-barrier topography in all simulations was fixed and specified by power laws with an exponent of 0.50 and constants of 0.15 and 0.10, respectively, in the transition from 0 m, MSL to a nominal island elevation of 1.5 m, MSL.

The dune was specified within the subaerial region spanning $x_d$ in Fig. 2. While the initial dune volume was held constant in all simulations, the dune profile takes on different shapes. The dune profile shape is based on the normal probability density function, so that different variances result in different shapes: a large variance corresponds to wider and shorter dunes (light orange in Fig. 2), while a small variance corresponds to narrower and taller dunes (dark orange in Fig. 2). The initial subaerial sediment volume within the region spanning $x_d$ is 406 m³/m (59% of total subaerial island volume $V$). Dune elevation, $z_d$, varied from 2.5 to 5.5 m, with lower heights corresponding to wider dunes and higher heights corresponding to narrower dunes. Finally, distance from the ocean shoreline to the initial profile’s dune crest, $x_c$, was varied from 159 to 249 m (Table 1).

These selected idealized profile conditions fall within the range of observed topographies along the United States Atlantic and Gulf or Mexico coasts (e.g., Cheng 2015; Cañizares and Irish 2008; USACE et al. 2022; NOAA 2022).

The horizontal computational grid resolution varies from 2 to 45 m, where finer resolution is used in areas with larger topographic gradients.

Sediment grain size was represented by a distribution characterized by the 50th and 90th percentile grain diameters, respectively $d_{50}$ and $d_{90}$ (Table 1). To understand the relative importance of grain size, we considered three different sands varying from fine ( $d_{50}=0.2$ mm) to coarse ( $d_{50}=0.6$ mm), which are characteristic of sandy barrier islands in the United States (e.g., Bocamazo et al. 2011; Feagin and Williams 2008). In each simulation, the grain size distribution was held uniform across the entire profile.

We assumed that belowground biomass acts to slow sediment mobility and, as such, can be simulated within XBeach by reducing equilibrium sediment concentration, $C_{\mathrm{eq}}$, when biomass is present. We make use of preexisting XBeach model parameters within XBeach’s governing equations; no new parameters or modules are introduced into the modeling framework. This is achieved within the model by introducing into the two-dimensional advection-diffusion equation an equilibrium sediment concentration factor $\theta$ (XBeach’s $sedcal$ parameter):where $C$ is the sediment concentration, $u$ and $v$ are the horizontal flow velocities in the $x$ and $y$ directions, respectively, $h$ is the water depth, and $D_h$ is the diffusivity. The factor $\theta$ varies from zero to one, where $\theta =1$ represents sediment alone (no biomass). When $\theta <1$, proportionately less sediment can be picked up and, thus, transported.

In our simulations, we made use of ten vertical sediment layers to approximate the impact of belowground biomass overlaying sediment. The top nine layers are each 5 cm deep and the tenth layer encompasses all sediment below the ninth layer. Sediment mobility was impeded in the top one or two sediment layers spanning $x_d$ in Fig. 2. In this parameter study, we assumed an impeded mobility layer thickness range of 5 to 10 cm, based on our own observations (Fig. 1) that the fibrous roots likely to provide the most resistance to erosion are more densely packed toward the surface (e.g., Schweiger and Schuettrumpf 2021). Field observations in New Jersey, USA (Charbonneau et al. 2016) and North Carolina, USA (Walker and Zinnert 2022), for example, showed that maximum root depth of the native American beach grass (Ammophila breviligulata) prevalent along the US Gulf of Mexico and Atlantic coasts extended no more than approximately 40 cm below the surface. Charbonneau et al.’s field photographs demonstrate a sparsening of fibrous roots with depth from the surface and with lateral distance from the aboveground shoots. However, others reported deeper Ammophila breviligulata maximum root depth (e.g., USDA 2013), and Charbonneau et al. and others have noted that nonnative grass species can have much deeper root depths. The impeded mobility layer thickness ${\mathrm{\Delta }}_{gl}$ was held uniform across the profile but varied from 5 to 10 cm, depending on the simulation (Table 1). In the initial profile’s impeded mobility layer in the biomass region $x_d$, $\theta$ was varied from 0.01 to 1.0 (Table 1). Beneath the impeded mobility layer, $\theta$ was initially fixed at 1 (no reduction in sediment mobility) to represent sand. Outside of the biomass region, $\theta$ was likewise initially fixed at 1. The reduced- $\theta$ material initially in the impeded mobility layer is free to move between vertical layers during the simulation, as this material is eroded and deposited, or covered up. In this regard, our approach is similar to that of van der Lugt et al. (2019) in that both approaches facilitate time-varying sediment mobility as belowground biomass is eroded or buried. Fig. S1 presents representative initial, interim, and final profiles depicting the redistribution over time of material initially in the top vertical layers.

While cohesion in the belowground biomass will certainly impact sediment slumping by avalanching, in our approach, we do not make any specific adjustments for this process. Here, we used XBeach critical wet and dry avalanching slopes within the standard range given in the XBeach manual, $\mathrm{wetslp}=0.1$ (similar to XBeach BOI value of 0.15, Rijkswaterstaat 2022) and $\mathrm{dryslp}=1.0$ (XBeach manual’s listed default). XBeach model parameter settings related to avalanching and sediment layers are provided in Table S1.

In using the combined effect of the model parameters $\theta$ and ${\mathrm{\Delta }}_{gl}$ on sediment mobility as an analogue for the impact of belowground biomass, we do not claim that the model setup exactly mimics real belowground biomass properties. From this perspective, the model parameters $\theta$ and ${\mathrm{\Delta }}_{gl}$ represent belowground biomass calibration parameters.

To assist in isolating the impact of dune shape and belowground biomass, idealized storm conditions were used in all simulations (Cheng 2015). Storm surge and wave conditions were selected to emphasize simulation of the Sallenger (2000) overwash and inundation regimes. This is a purely hypothetical storm that simplifies both surge height and wave height variation with time. Real storms may have shorter or longer durations and may exhibit some skewness around the surge peak, for example. In all simulations, storm significant wave height and peak period at the offshore boundary (25-m depth) are respectively 8 m and 10 s, where a unidirectional onshore JONSWAP spectrum is employed. Idealized storm surge was specified using the positive portion of a sine curve over a 6-hour duration, starting from 1.0 m, MSL then increasing to a peak of either 2.5 or 3.0 m, MSL before decreasing to 1.0 m, MSL (Table 1).

Simulations were planned for all unique input parameter combinations in Table 1 for a total of 6,480 attempted simulations. Of these, 36 simulations (0.6%) failed leaving 6,444 completed simulations for analysis. For each completed simulation, the final bed elevation profile was saved for analysis. From these profiles, bed elevation change with respect to the initial profile is determined then used to quantify the following.

These output parameters were examined to determine their sensitivity and trends with respect to model input parameters, with emphasis on the influence of impeded sediment mobility. We used partial correlation analysis to isolate the influence of one input parameter while controlling for other input parameters. We further selected Spearman correlation to quantify correlation while allowing for nonlinear monotonic trends. For discussion, we interpret Spearman partial correlation coefficients $r$ as follows.

Fig. 3 depicts simulated bed elevation change along the cross-shore profile as a function of equilibrium sediment concentration factor $\theta$ for selected dune positions ( $x_c$) and elevations $z_d$ when peak surge is 2.5 m, MSL, median grain diameter is $d_{50}=0.4$ mm, and impeded mobility layer thickness is ${\mathrm{\Delta }}_{gl}=10$ cm. Fig. 4 likewise depicts the same information for the 3.0-m, MSL peak surge case. Regardless of dune position and elevation, lower $\theta$ results in less sediment movement along the profile, both in terms of magnitude of vertical change and lateral extent of bed change. With the exception of the highest dune elevation case ( $z_d=5.5$ m, MSL) during the lower-surge event, for a given dune elevation, the amount of erosion at the dune crest is similar for the narrower [panels (a, c, e, g, and i)] and wider [panels (b, d, f, h, and j)] distance to dune crest cases. Yet, lateral extent of overwash deposit shifts landward (to the right) of the order of the difference in distance to the dune crest. The dune elevation of $z_d=5.5$ m, MSL during the lower-surge event varies from predominantly overwash regime to predominantly inundation regime, depending on the factor $\theta$ and distance to dune crest $x_c$, with the wider $x_c$ case exhibiting more limited dune elevation erosion.

Respectively for the 2.5- and 3.0-m, MSL peak surge cases, Figs. 5 and 6 depict simulated bed elevation change along the cross-shore profile as a function of equilibrium sediment concentration factor $\theta$ for each of the three simulated grain size distributions and each of the two simulated impeded mobility layer thicknesses when the initial dune elevation is $z_d=4.2$ m, MSL and distance to the dune crest is $x_c=159$ m. As anticipated, sediment movement increases with decreasing sand grain diameter. However, a doubling of the impeded sediment layer thickness from 5 to 10 cm has a less dramatic influence on sediment movement; its influence is most pronounced for very low $\theta$.

Table 2 lists Spearman partial correlation coefficients $r$ between input and output parameters when all 6,444 simulations are considered. All input–output parameter pairs result in statistically significant correlation, with $p$-values less than 0.05. As anticipated, all output parameters exhibit moderate to strong correlation with peak surge, with positive correlation with shoreline change $\mathrm{\Delta }{x}_o$, total subaerial eroded volume $V$, and overwash extent $x_{\mathrm{ow}}$ and negative correlation with final dune elevation $z_{d\text{-final}}$.

With respect to the two topographic characteristics, all output parameters are moderately correlated with initial dune elevation $z_d$. Because initial subaerial profile volume is held constant in all simulations, a positive correlation with initial dune elevation also indicates a negative correlation with dune width, and vice versa. In contrast, all but shoreline change $\mathrm{\Delta }{x}_o$ are uncorrelated (defined as $|r|<0.25$) with distance to dune crest $x_c$. Shoreline change exhibits moderate positive correlation with $x_c$.

With respect to sediment and sediment mobility characteristics, output parameters are moderately to strongly correlated with the factor $\theta$ and with sand grain size $d_{50}$. In contrast, output parameters are uncorrelated to weakly correlated with impeded mobility layer thickness ${\mathrm{\Delta }}_{gl}$.

Overall, this Spearman partial correlation analysis on the full 6,444 simulation set suggests that sediment movement under the overwash and inundation regimes is dominated primarily by sediment grain size, the factor $\theta$, and initial dune elevation. While the presence of the impeded mobility layer has a significant impact on sediment movement, as indicated by strong correlation with $\theta$, the thickness of the impeded mobility layer ${\mathrm{\Delta }}_{gl}$ only marginally influences sediment movement. The results further suggest that for these two storm regimes, distance to the dune crest does not measurably impact sediment movement.

Next, we investigate the range of output parameter responses as a function of the commonly accepted primary drivers of sediment change during inundation and overwash, initial dune height, and peak surge, while considering all other input parameters as sources of uncertainty. Fig. 7 depicts box plots of the output parameters when the simulations are grouped by like initial dune elevation (642 to 648 simulations per box). Total subaerial volume eroded $V$ [Fig. 7(e and f)] exhibits the most variation across a range of initial dune elevations, often varying by an order of magnitude. As anticipated, median $V$ decreases with increasing initial dune height, exhibiting a stronger trend with initial dune height for the lower-surge case [Fig. 7(e)].

Shoreline change [Fig. 7(c and d)] likewise exhibits variation across a range of initial dune elevations though the magnitude of change of the order of 20 m or less is relatively small with respect to distance to the dune. Notably, the direction of shoreline change generally shifts to accretion (negative) for the lower surge case [Fig. 7(c)] when dune elevation is 4.7 m or higher, corresponding to cases predominantly in the overwash regime.

Overwash extent [Fig. 7(g and h)] tends to be fairly consistent across simulations of the same initial dune height, except in the 4.7-m dune elevation case for the lower surge, when simulations span from being overwash-dominated to inundation-dominated. Final dune elevation [Fig. 7(a and b)] is likewise constrained to a small range close to the base island height of 1.5 m in all but the two highest initial dune elevation cases for the lower surge [Fig. 7(g)]. As with overwash extent, this wider range in final dune elevation at the higher initial dune heights reflects simulations spanning from overwash-dominated to inundation-dominated.

Fig. 8 depicts Spearman partial correlation between input and output parameters when data are grouped by like initial dune elevation and surge magnitude. With respect to sediment mobility, regardless of initial dune height, the total subaerial eroded volume $V$ [Fig. 8(e and f)] exhibits very strong correlation with both sediment grain size $d_{50}$ and the factor $\theta$, and is moderately correlated with impeded mobility layer thickness ${\mathrm{\Delta }}_{gl}$. Except for the lower surge case [Fig. 8(a, c, e, and g)] when dune height is 4.7 m, MSL or higher, final dune elevation $z_{d\text{-final}}$ [Fig. 8(a and b)] and overwash extent [Fig. 7(g and h)] likewise are strongly correlated with $d_{50}$ and $\theta$; these output parameters are moderately to weakly correlated with ${\mathrm{\Delta }}_{gl}$.

Consistently, across initial dune heights, shoreline change $\mathrm{\Delta }{x}_o$ exhibits strongest correlation with grain size $d_{50}$ [Fig. 8(c and d)]. Shoreline change mostly exhibits moderate positive correlation with the factor $\theta$. Yet, for the highest dune and lowest surge case, shoreline change’s correlation with $\theta$ becomes weakly negative. Of the four output parameters, only shoreline change exhibits consistent dependence on distance to the dune crest $x_c$; here, correlation is moderate to strong for most initial dune height and storm condition combinations. Overwash extent’s dependence on $x_c$ varies from no correlation to highly correlated, with the highest correlation corresponding to the overwash-dominated case with low surge and high dune height. Final dune elevation and total subaerial eroded volume both exhibit weak to no correlation with $x_c$.

The following sections discuss the nature of the relationship between the four analyzed output parameters and correlated input parameters.

Fig. 9 depicts total subaerial eroded volume for the 2.5-m surge case when distance to the dune crest is $x_c=159$ m, as a function of the factor $\theta$ and median grain diameter $d_{50}$; results are shown for different initial dune elevations $z_d$ and for both impeded mobility layer thicknesses. Fig. 10 provides the same information when peak surge is 3.0 m. In general, higher and narrower dunes (higher $z_d$) result in less total erosion, although the initial subaerial island volume is kept constant.

Regardless of storm surge intensity and initial dune elevation, a lower $\theta$ results in less volume of total erosion. Decreasing sediment mobility has a greater impact on total volume eroded when sediment size is small; volume eroded is more variable with $\theta$ when median diameter is smaller and less variable when diameter is larger. In other words, the largest gradient in total volume generally occurs when both $\theta$ and $d_{50}$ are small while the smallest gradients are associated with higher $\theta$ and larger $d_{50}$. The magnitude of these gradients increase with decreasing initial dune elevation, indicating that the impact of impeded mobility is more significant when the dune is narrower and wider, or alternatively, is more deeply inundated. The pattern of total eroded volume in the $\theta$- $d_{50}$ space for the lower surge and highest dune elevation case [Figs. 9(a and b)] differs from the other surge and dune elevation cases, where in most simulations, the higher dune sufficiently protected against full inundation (i.e., simulations were overwash-dominated). Here, volume eroded approaches a linear trend with sediment size.

Finally, as indicated by the correlation analysis, total eroded volume is marginally impacted by impeded mobility layer thickness ${\mathrm{\Delta }}_{gl}$. Comparing panels (a, c, e, g, i) ( ${\mathrm{\Delta }}_{gl}=5$ cm) and panels (b, d, f, h, j) ( ${\mathrm{\Delta }}_{gl}=10$ cm) in Figs. 9 and 10 reveals that the relationship between eroded volume and $\theta$ and $d_{50}$ is very similar.

The influence of sediment mobility and sediment size on final dune elevation $z_{d\text{-final}}$ and on overwash extent $x_{\mathrm{ow}}$ are similar, as in Figs. 11 and 12 for simulations with initial dune height $z_d=5.5$ m, MSL and distance to dune crest $x_c=169$ m. This similarity is expected, given that more dune erosion means more sediment is available to be moved into the back-barrier region. Our simulation results yield a Spearman correlation coefficient of $r=-0.93$ between these two output parameters.

Simulations depicted in Figs. 11 and 12 range from those that are predominantly overwash-regime dominated [panels (a and b)] to those that are predominantly inundation-regime dominated [panels (c and d)]. Panels (c and d) in these figures depict that when the island is unindated, the influence of the factor $\theta$ on $z_{d\text{-final}}$ and $x_{\mathrm{ow}}$ is limited to very small $\theta$, of the order of 0.1 or smaller. For larger $\theta$, $z_{d\text{-final}}$ and $x_{\mathrm{ow}}$ are nearly constant for a given sediment grain diameter, and variation with grain diameter is limited to a relatively small range.

In contrast, Figs. 11(a and b) exhibit strong variation in $z_{d\text{-final}}$ across all values of $\theta$. The $z_{d\text{-final}}$ is most sensitive to $d_{50}$ when $\theta$ is larger, being close to constant when $\theta$ is less than approximately 0.2. The 4.0–4.3 m, MSL contour interval represents the transition from overwash-dominated to inundation-dominated conditions, exhibiting a much stronger dependence on grain diameter as the dune lowers below 4.3 m, MSL. This threshold flow condition is also apparent in Figs. 12(a and b), where overwash extent likewise rapidly increases with decreasing $d_{50}$ and increasing $\theta$. Prior to reaching this threshold, overwash extent is limited to the vicinity of the dune.

While the shoreline’s position is seaward of the belowground biomass region, shoreline change exhibits sensitivity to the presence of the impeded mobility layer, as seen in Fig. 13 when the initial dune elevation is $z_d=5.5$ m and impeded mobility layer thickness is ${\mathrm{\Delta }}_{gl}=10$ cm. Simulations shown for the lower surge case in Figs. 13(a and b) are generally overwash-regime dominant. The resulting shoreline response is accretion, as a result of mobilized sediment being redeposited seaward of the dune. In contrast, for inundation-dominant simulations, as is the case for the high surge [panels (c and d)], the resulting shoreline response is erosion and indicates a net loss of sediment from the region seaward of the dune.

Decreasing $\theta$ in the biomass region leads to less shoreline erosion during inundation events and less accretion during overwash events – reflecting an overall reduction in mobile sediment with lower $\theta$. While increasing sediment grain size $d_{50}$ likewise leads to less shoreline erosion during inundation events, increasing grain size leads to more shoreline accretion during overwash events.

Of the output parameters we considered, shoreline change is most sensitive to initial dune position. As distance from the shoreline to the dune crest $x_c$ increases from narrower [panels (a and c)] to wider [panels (b and d)], less sediment is deposited at the shoreline during overwash events while more sediment is removed from the shoreline during inundation events.

While our results reinforce that initial dune elevation, hydrodynamic forcing, and sediment grain diameter remain dominant factors driving barrier-island change during storm overwash and inundation, our results also demonstrate that impedance of sediment mobility, as is the case in the presence of belowground biomass, can play a significant role in barrier-island change. Over the parameter ranges considered herein, the factor $\theta$ is potentially as important an indicator of eroded volume, dune elevation change, and overwash extent as is initial dune elevation and grain diameter (Table 2). The practical engineering significance of this result lies in quantitatively demonstrating the importance of designing and maintaining vegetation elements in conjunction with traditional beach nourishment topographic parameters.

As observed by van der Lugt et al. (2019) and Schweiger and Schuettrumpf (2021), our findings demonstrate the significance of belowground biomass in reducing the severity of erosion and its potentially important role in the transition between overwash and inundation regimes. Though our approach differs from van der Lugt et al.’s (2019) in that we impede sediment mobility rather than reduce flow velocity through roughness, both approaches allow for time-dependent variation in sediment erosion and deposition as a result of differing bed characteristics. While dune height remains a dominant driver of when and whether the flow regime transitions from overwash to inundation, this bed-characteristic time dependence is important. When total water level coincides with dynamic dune elevation, belowground biomass may well become the dominant factor in determining whether flow remains within the overwash regime or transitions to the inundation regime, particularly when threshold conditions occur near the water-level peak.

This study represents the first attempt to quantify barrier-island morphological response over a range of impeded mobility, grain diameter, and dune elevation. Our results show that the interplay between impeded sediment mobility and grain diameter is complex (Figs. 9–13). Yet, in nature, both dune vegetation and sediment characteristics vary not only at regional scales but also at local scales of the order of our finest model grid resolution (2 m) and less. Vegetation in particular varies over time as well. Such local-scale observational data are not readily obtained or available, though initiatives like the citizen-science sand grain project SandSnap (e.g., McFall et al. 2022) show promise for the future. In the meantime, our results may best be employed to quantify aggregate uncertainty in expected morphological response and in ensemble simulation applications (e.g., Monte Carlo) to characterize the range and likelihood of various morphological responses.

As is introduced in the Methods section, the belowground biomass-induced mobility reduction factor $\theta$ is a scaling factor applied to the equilibrium sediment concentration. Lower $\theta$ values reduce the equilibrium concentration within the top layers of bed sediment of specified thickness ${\mathrm{\Delta }}_{gl}$, and thereby limit the amount of sediment picked up by the flow. Our results show a strong dependence on $\theta$, but a relatively weak dependence on ${\mathrm{\Delta }}_{gl}$. This relatively weak influence of changes in ${\mathrm{\Delta }}_{gl}$ on morphological response is surprising. To further explore our finding regarding the relatively low importance of the impeded mobility layer thickness ${\mathrm{\Delta }}_{gl}$, we conducted a computational sensitivity analysis for the $z_d=4.2$ and $x_d=159$ m dune case when $d_{50}=0.4$ mm and $\theta =0.05$, 0.1, 0.2, 0.3, and 0.5. We varied ${\mathrm{\Delta }}_{gl}$ at 5-cm increments from 5 to 45 cm, with the deepest case being in the range of van der Lugt et al.’s (2019) somewhat arbitrarily selected belowground biomass XBeach calibration depth. Results for these 45 simulations are in Fig. 14. Of the four outputs, total subaerial eroded volume is the most sensitive to ${\mathrm{\Delta }}_{gl}$, where its influence is most pronounced for $\theta \le 0.1$. For $\theta \ge 0.2$, the interquartile range of eroded volume is relatively small, within $\pm 13$ $m^3/m$ (13%) of the median result. Likewise, final dune elevation $z_{d\text{-final}}$ and overwash extent $x_{\mathrm{ow}}$ are only impacted when $\theta \le 0.1$. When $\theta \ge 0.2$, respectively for $z_{d\text{-final}}$ and $x_{\mathrm{ow}}$, the interquartile ranges are within $\pm 0.10$ and $\pm 10$ m of the median result; these differences are of the order of the numerical error and thus can be considered negligible. Variability in shoreline change $\mathrm{\Delta }{x}_o$ across all values of $\theta$ is likewise negligible, with changes of no more than 5 m falling within the range of numerical error. Under the assumption that fibrous root density thins with depth below the surface, this sensitivity analysis largely supports our earlier finding that results are relatively insensitive to impeded mobility layer thickness ${\mathrm{\Delta }}_{gl}$. However, in cases where very dense, fibrous roots extend deep into the dune, impeded mobility layer thickness could become important.

It should be noted that $\theta$ and ${\mathrm{\Delta }}_{gl}$ are not the only parameters in XBeach that affect vegetated flow and sediment transport. Instead, there is a whole group of parameters to define very detailed vegetation characteristics that influence vegetated flow, such as, vegetation stem diameter, vegetation density, and vegetation height per vegetation layer, to name a few. Using these parameters will lead to damping of waves which, in turn, will affect sediment transport (Van Rooijen et al. 2015).

The dissipation of waves mentioned above only applies to stationary vegetation and does not account for any vegetated sediments entrained by the flow. In comparison, a lower-than-one value for the factor $\theta$ not only affects erosion, but it will accelerate the deposition of sediment in the flow as well. In reality, the effect of vegetation on sediment transport is much more complicated and our knowledge is currently too limited to be able to fully model vegetated sediment transport as a process. For instance, the shape of vegetation itself – both its belowground root system and aboveground stem and leaf system – is a challenge for transport modeling; vegetation could capture flowing sediments, and vegetation could greatly alter the properties of surrounding sediments by making them cohesive and increasing the effective grain size (Feagin et al. 2019).

Practical application of the model parameters $\theta$ and ${\mathrm{\Delta }}_{gl}$ to model specific belowground biomass necessitates calibration to field and laboratory data. There is currently no method for direct mapping of certain $\theta$ and ${\mathrm{\Delta }}_{gl}$ values to a specific combination of vegetation species, maturity, and health. However, such relationships can be established, given sufficient field or experimental data of vegetation characteristics, flow conditions, and sediment transport. van der Lugt et al. (2019) demonstrated successful regional-scale calibration of their dynamic bed roughness model using initial roughness specified using spatially variable land cover classification and using field observations of overwash and breaching during hurricanes in Florida and New York, USA. A calibration exercise similar to that undertaken by van der Lugt et al.’s (2019) is a logical next step in mapping the parameter $\theta$ to land cover classification for practical application.

As performed by Schweiger and Schuettrumpf (2021), controlled laboratory experiments may likewise be used to constrain or calibrate model parameters. For example, as shown in Fig. 15, the recent work of Bryant et al. (2019) reveals that our simulations fall within the range of their belowground biomass laboratory results. The gray shaded area in the figure spans the range of Bryant et al. (2019) dune erosion results for the case where artificial belowground biomass is introduced uniformly throughout the dune and the dune overwashes. Our simulations suggest that setting the factor $\theta$ between 0.1 and 0.3 would lead to reasonable XBeach results for this experimental case.

In summary, our use of $\theta$ and ${\mathrm{\Delta }}_{gl}$ to represent belowground biomass is an intuitive approach that works within XBeach’s existing governing equations and settings to approximate the impact of dense, fibrous roots on sediment mobility (see Methods section). Yet, this is a rather tentative approach to simulating the influence of vegetation on sediment transport. As such, this approach is currently difficult to verify and validate. Nonetheless, its simplicity allows us to represent the influence of belowground biomass in an intuitive way with just two model parameters, $\theta$ and ${\mathrm{\Delta }}_{gl}$, which makes possible a general parameter study that involves many other variables.

In this study, we construct different dune shapes based on the probability density function of normal distributions. This approach allows the change of dune shape (height and width) while retaining a constant volume within the dune. However, a drawback of this approach is that for larger initial dune elevations ( $z_d$), the actual dune width is smaller than the defined belowground biomass region width $x_d$. For example, the transition from the berm to the dune (dune toe) is typically determined by slope and would thus result in a narrower dune width. Our assumption that the biomass region $x_b$ is constant regardless of $z_d$ leads to impeded sediment motion outside of the actual dune width when $z_d$ is large.

This study is the first parameter study aimed at investigating the impact of belowground vegetation biomass on barrier island response during storms. We have considered parameters including dune elevation, distance to the dune crest, sediment grain size, and belowground biomass-induced sediment mobility reduction. Our simulation results have demonstrated how each individual parameter affects subaerial erosion and deposition, and more importantly, revealed the complex interplay between parameters. Our results can also be employed in barrier island storm hazard prediction and protection efforts. Furthermore, our approach of implementing the influence of belowground biomass using existing XBeach model settings and governing equations may be readily adopted by practitioners, as it negates the need to implement new modules that are potentially more computationally intensive. Additional insights can be gained by expanding the number of input parameters (six) considered along with expanding the number of distinct values considered for each parameter. Such an expanded parameter study needs to take into account variations in, for instance, horizontal and vertical distribution of belowground biomass and aboveground stems and leaves, dune volume, elevation of berm and back barrier, and storm conditions (wave height, surge magnitude, and duration) – to include inducing conditions leading to Sallenger’s (2000) collision regime. Again, the purpose of this study is not focused on reproducing any specific storm case at any location, but to provide coastal engineers and scientists with a tool to better understand the sediment transport on vegetated barrier islands during storms. For dune design and restoration efforts, our results help point out the locations where most work is needed, and to what extent, if vegetation is involved in the dune design. This study has laid a solid groundwork for future studies, with the ultimate goal of improving our understanding of barrier island dynamics and better protecting the coastal communities from storms.