Considering the Effect of Land-Based Biomass on Dune Erosion Volumes in Large-Scale Numerical Modeling

This paper presents and validates a novel root model which accounts for the effect of belowground biomass on dune erosion volumes in XBeach, based on a small-scale wave flume experiment that was translated to a larger scale. A 1D-XBeach model was calibrated by using control runs considering a dune without vegetation. Despite calibration, a general model–data mismatch was observed in terms of overestimated erosion volumes around the waterline. Furthermore, the prediction of overwash had to be induced by increasing the maximum nearshore wave height within the XBeach simulation. Subsequently, applying the root model resulted in a good agreement with the belowground biomass cases, and the consideration of spatially varying rooting depths further improved the results. Predictions of the root model while using locally increased friction coefficients were in line with the aboveground and belowground biomass cases. However, the effect of the root model on the erosion predictions varied among the hydrodynamic conditions, so further improvements are required. Therefore, future research should focus on quantifying the effects of land-based biomass and individual plant characteristics, such as root density, on dune erodibility at large scales, along with their influences on the temporal evolution of dune scarping and avalanching.


Introduction
Coastal dunes are distributed worldwide and acknowledged for their wide range of functions, including their contributions to coastal defense, ecological diversity and socioeconomic services, such as recreation and tourism [1,2]. Concerning coastal protection, dunes serve as a natural barrier at the boundary between land and sea. Furthermore, they have been recommended as one example of a more sustainable, cost-effective and ecologically sound coastal protection measure than conventional hard structures [3][4][5]. In this regard, dune vegetation was found to have an especially beneficial effect on the resilience of dunes against storm-induced erosion [6][7][8][9][10]. Due to the globally observed climaticallydriven increase in dune vegetation in the period between 1984 and 2017 [11] and the expected increasing frequency of sea level extremes until the end of this century [12][13][14][15][16], detailed knowledge about wave-driven erosion of vegetated beach-dune systems is important with respect to global climate change. Owing to its wide use in coastal management, understanding the influence of vegetation on nearshore hydrodynamic and morphodynamic processes is also essential concerning coastal numerical modeling.
The process-based XBeach model [17] has been proven to accurately predict storminduced erosion and dune breaching [18][19][20][21][22][23]. With respect to vegetation, the model was extended to account for the damping of short waves and infragravity waves and mean flow [24,25], and the interaction between waves and vegetation [26]. From a morphodynamic perspective, the effect of vegetation is often considered by locally increasing the bed friction coefficient [22,[27][28][29], which was improved by including a dynamic response of local bed friction coefficients in vegetated areas [30]. In addition, Bendoni et al. [31] improved the XBeach model in order to model the erosive processes of saltmarsh vegetation.

XBeach
XBeach [17] is a two-dimensional (2DH) numerical model which solves equations for flow, wave motion, sediment transport and bed evolution. While short-wave motion is obtained from a time-dependent version of the wave-action balance equation on the wave group time scale (surf beat mode), low-frequency and mean flows are calculated based on the depth-averaged shallow-water equations. These are cast into a depth-averaged generalized Lagrangian mean (GLM) formulation to account for wave-induced mass flux and subsequent return flow.
The bed friction associated with low-frequency and mean flow is considered by the bed shear stress (τ) approach according to Ruessink et al. [37]: with ρ being the density of water, u E and v E the Eulerian velocities and u rms the near-bed short-wave orbital velocity. The dimensional friction coefficient c f can be determined, among others, by the Manning coefficient n: with g being the acceleration due to gravity and h the water depth. In order to account for spatial and temporal variations, a dynamic roughness module [30] varies the Manning roughness coefficients in vegetated areas according to: n =    n sand + n veg − n sand ·min max 0, d * +∆zb d * , 1 , −d * ≤ ∆zb < 0 n sand + n veg − n sand ·min max 0, h * +∆zb h * , 1 , ∆zb > 0 3 of 21 with n sand/veg being the Manning coefficients associated with sand and vegetation, d* and h* the critical depths for erosion and deposition and ∆zb the cumulative deposition/erosion which is negative in the case of erosion. Sediment transport is modeled with the aid of a depth-averaged advection-diffusion equation [38], where the actual sediment concentration C is calculated by the mismatch with an equilibrium sediment concentration C eq . The latter can be calculated by multiple sediment transport formulations, of which one example is the Soulsby-Van Rijn equations [39,40], in which the equilibrium sediment concentration for bed and suspended load is calculated according to: with A sb and A ss being the bed load and suspended load coefficients, v mg the Eulerian velocity magnitude, u rms,2 the adjusted near-bed short-wave orbital velocity for wave breaking induced turbulence, C d a drag coefficient and U cr the critical velocity that defines at which depth-averaged velocity sediment motion is initiated. The actual sediment concentration is used for the calculation of sediment transport rates, which in turn serve as input for the calculation of the bed-level update. Finally, an avalanching mechanism accounts for the slumping of sandy material, which exchanges sediment between two adjacent cells as long as a critical slope between these cells is exceeded. By distinguishing between wet and dry cells, the avalanching mechanism considers that wet cells are more prone to slumping. A full description of the XBeach model can be found in Roelvink et al. [17] and in the online manual (https://xbeach.readthedocs.io/en/latest/xbeach_manual.html#id83, accessed on 6 July 2021).

Root Model
The fundamental idea of the root model is based on the effect of fibrous roots on sandy soils by increasing the resistance of the top (rooted) soil against concentrated flow erosion due to additional root cohesion [41]. Figure 1 shows the basic principle of the root model. Within the morphevolution module of the XBeach model, a root cohesion term increases the critical erosion velocity where belowground biomass is present until the cumulative erosion (∆zb) exceeds a user-defined rooting depth (z root ) ( Figure 1a): , z root ≤ ∆zb < 0 (5) with ρ being the density of water. The user-defined root cohesion C r (kN/m 2 ) is the product of the root tensile strength t R and the root area ratio RAR, which are both plant-specific [42]. An additional root cohesion coefficient (rcc) is used to steer the model, which can either be constant or linearly vary with the erosion depth ( Figure 1b): , 0 , 1 ·rcc 0 , dynamic mode (6) with rcc 0 being the initial root cohesion coefficient which is zero in areas without belowground biomass. In comparison to its first version, the root model was extended to account for spatially varying rooting depths (z root (x)), which are provided through an external file that has the same format as the hardlayer file (Figure 1c). Similarly to its application to smaller scales [32], the root cohesion term was normalized to one so that applying rcc = 1 corresponds to an increase of U cr by U cr,root = 1 m/s. This was done to simplify the calibra-tion process in this study and due to the use of substitutes (coir fibers) with unknown t R and RAR instead of natural vegetation.
ci. Eng. 2021, 9, x FOR PEER REVIEW 4 of 21 external file that has the same format as the hardlayer file (Figure 1c). Similarly to its application to smaller scales [32], the root cohesion term was normalized to one so that applying rcc = 1 corresponds to an increase of Ucr by Ucr,root = 1 m/s. This was done to simplify the calibration process in this study and due to the use of substitutes (coir fibers) with unknown tR and RAR instead of natural vegetation. For a detailed description of the root model, the reader is referred to Schweiger and Schuettrumpf [32].  simulated by the root model until the cumulated erosion (Δzb) exceeds a user-defined rooting depth (zroot). The root cohesion term is steered by a root calibration coefficient (rcc) which can be constant or linearly decreasing with the erosion depth (dynamic mode). (c) In addition to the first version [32], it is now possible to use spatially varying rooting depths.

Model Scaling
The aim of this study was to validate the root model on a larger scale problem by upscaling the small-scale model setup of Bryant et al. [34]. In general, scaling of laboratory experiments requires the correct representation of the physical processes in nature so that dimensionless numbers characterizing these processes (e.g., Froude number, Reynolds number and Irribaren number) are the same for the prototype and model [43]. The most common approach is the use of scaling laws, where the ratio between prototype and model is defined by the scale parameter n = pp/pm. Here, pp and pm are the parameter values in prototype and in the (laboratory) model [43]. In this regard, hydrodynamic parameters are generally scaled according to Froude [43][44][45][46][47]:  For a detailed description of the root model, the reader is referred to Schweiger and Schuettrumpf [32].

Model Scaling
The aim of this study was to validate the root model on a larger scale problem by upscaling the small-scale model setup of Bryant et al. [34]. In general, scaling of laboratory experiments requires the correct representation of the physical processes in nature so that dimensionless numbers characterizing these processes (e.g., Froude number, Reynolds number and Irribaren number) are the same for the prototype and model [43]. The most common approach is the use of scaling laws, where the ratio between prototype and model is defined by the scale parameter n = p p /p m . Here, p p and p m are the parameter values in prototype and in the (laboratory) model [43]. In this regard, hydrodynamic parameters are generally scaled according to Froude [43][44][45][46][47]: where H is the wave height, L is the wave length, h is the water depth, T is the wave period, t is the simulation time and u is the wave orbital velocity. Concerning physical geometry, the distortion scale ratio is used [43]: with n g being the distortion scale, n l the length ratio and n h the height ratio between prototype and model. In coastal modeling, these are usually represented by the length and the height of the beach profile being analyzed [47]. From n g = 1, it follows that the physical geometry of the (beach-dune) model is undistorted. If the same applies to Froude scaling, two models will automatically have the same Irribaren number and wave steepness (H/L) [47]. With respect to sediment characteristics, however, additional dimensionless parameters need to be considered. Examples are the dimensionless fall velocity (Dean number), the Shield's number and the Reynolds number [43,47]. From these follows the scaling parameter n D50 , for which various scaling laws exist, as elaborated in van Rijn et al. [43]. Concerning coastal dune erosion, the following scaling relation is proposed: with D 50 being the medium sediment diameter and (s − 1) the relative density. In case of sand and an undistorted scale, Equation (9) reduces to n D50 = n h 0.56 . However, as discussed in Bayle et al. [47], sand grain size scaling always results in scaling errors: If, for instance, Froude scaling is maintained and sand is correctly scaled, then the Shield's number and the Dean number are maintained as well, whereas the grain Reynolds number is not. As a result, sand scaling is often not performed due to the permanent availability of the same type of sand and potential cohesive properties for small-scale models.

Model Setup
A good example of XBeach large-scale modeling is the application to Deltaflume test T04, which was also used for initial calibration of the XBeach model (see Roelvink et al. [17] for the details). A small dune with a height of approx. 1.65 m was located in front of a larger dune and eroded due to collision and overwash [35]. Within the scope of this study, it was therefore chosen to upscale the small-scale 1D-model used in Schweiger and Schuettrumpf [32] to match the dune dimensions of Deltaflume test T04. From the initial small-scale dune height of h D = 0.223 m above still water level follows an increase in height by approximately a factor of SF~15/2 = 7.5 = n h . Applying an equal length ratio (n l = 7.5) results in a geometrical undistorted model (n g = 1).
The upscaled 1D-model with a length of approximately 475 m is shown in Figure 2. The origin of the local coordinate system equals the starting point of profile measurements (x = 0 m) and the initial lower water level (zb = 0 m), which was 225 cm above the flume floor (see Table 1). To reduce the computational cost, a fast 1D (ny = 0) nonequidistant grid was generated with dx max = 10 m at the offshore boundary and a minimal grid size of dx min = 0.1 m in the beach-dune area, which was chosen based on preliminary testing (nx = 253). The model setup was similar to that in Schweiger and Schuettrumpf [32]. Among others, spatially varying friction coefficients (Manning) were applied with n = 0.03 s/m 1/3 in the sandy section of the flume. For the concrete section, a Manning coefficient of n = 0.01 s/m 1/3 was used and hard structures enabled. With respect to sediment characteristics, grain size scaling would result in a medium diameter of D 50 = 0.46 mm (Equation (9)). However, the initial value of D 50 = 0.15 mm was maintained due to the general sediment scaling issues mentioned before (see Section 2.3). The four different hydrodynamic boundary conditions (shallow collision (SC1), shallow overwash (SO), deep collision (DC) and deep overwash (DO)) were scaled according to Froude. While the constant water levels (tideloc = 0) and Hm0 wave heights were scaled by SF = 7.5, the peak wave period Tp and the duration of each wave burst were scaled according to SF −0.5 . The wave boundary conditions were applied as JONSWAP spectra (γ = 3.3) with rt (duration of each spectrum at the offshore boundary) equaling the duration of one wave burst. The resulting simulation times were tstop = 1,100 s for DO (only one wave burst), tstop = 3,300 s for SO and tstop = 9,900 s for the C-conditions. Similarly to the application to small-scale, the default spin-up time of wave boundary conditions was reduced from taper = 100 s to zero in order to prevent underestimation of the measured wave heights during the first wave burst [32]. An overview of the upscaled hydrodynamic boundary conditions is given in Table 1. Table 1. Upscaled (hydrodynamic) boundary conditions according to Froude with still water level above the flume ground (SWL), zero-moment wave height Hm0, peak wave period Tp and the number and durations of wave bursts. The latter two define the applied simulation times (tstop).

Hydrodynamic Calibration
The model was hydrodynamically calibrated (morphology = 0) by comparing the upscaled measured Hrms,hf wave heights of each wave burst with the time-averaged spatial output Hmean over one wave burst (tintm = rt). The XBeach default settings were applied with gamma = 0.5, and thus, the same settings as for the small-scale applications [32]. Varying the breaker index between 0.45 < γ < 0.7 did not improve the overall results, so that a value of gamma = 0.5 was used for the remainder of the investigation.   The four different hydrodynamic boundary conditions (shallow collision (SC1), shallow overwash (SO), deep collision (DC) and deep overwash (DO)) were scaled according to Froude. While the constant water levels (tideloc = 0) and H m0 wave heights were scaled by SF = 7.5, the peak wave period T p and the duration of each wave burst were scaled according to SF −0.5 . The wave boundary conditions were applied as JONSWAP spectra (γ = 3.3) with rt (duration of each spectrum at the offshore boundary) equaling the duration of one wave burst. The resulting simulation times were tstop = 1,100 s for DO (only one wave burst), tstop = 3,300 s for SO and tstop = 9,900 s for the C-conditions. Similarly to the application to small-scale, the default spin-up time of wave boundary conditions was reduced from taper = 100 s to zero in order to prevent underestimation of the measured wave heights during the first wave burst [32]. An overview of the upscaled hydrodynamic boundary conditions is given in Table 1.

Hydrodynamic Calibration
The model was hydrodynamically calibrated (morphology = 0) by comparing the upscaled measured H rms,hf wave heights of each wave burst with the time-averaged spatial output Hmean over one wave burst (tintm = rt). The XBeach default settings were applied with gamma = 0.5, and thus, the same settings as for the small-scale applications [32]. Varying the breaker index between 0.45 < γ < 0.7 did not improve the overall results, so that a value of gamma = 0.5 was used for the remainder of the investigation. Figure  The comparison with the upscaled measured Hrms,hf wave heights shows good overall agreement for each hydrodynamic condition, with the bias varying between −0.3 cm (DC) and 5.77 cm (SO). The most agreement with the observation was achieved for the overwash cases with R² = 0.78 for SO and R² = 0.82 for DO. With respect to nearshore wave heights, the model underestimated the maximum observed wave height, except for SO. Overall, the results are very well in line with the model's small-scale application [32], as similar statistics were derived (see Table 2), which proves that the translation to a larger scale is valid with respect to hydrodynamics.  Table 2. Comparison of wave statistics between the small-scale (adapted from [32]) and large-scale (this study) applications based on the coefficient of determination (R²) and the bias for each hydrodynamic condition.

Morphodynamic Calibration
The morphodynamic calibration was performed based on the upscaled control cases and by reusing the boundary conditions generated during the hydrodynamic calibration (wbctype = reuse). Due to the similar model performance in small-scale and large-scale situations in terms of hydrodynamics, the final morphodynamic parameter settings of Schweiger and Schuettrumpf [32] were used as initial settings in this study; and the critical avalanching slopes were set to wetslp = 0.2 and dryslp = 0.7. However, the depthscale parameter was set to SF = 2 to be consistent with the upscaled model setup (1:2). As a result, the effected morphodynamic parameters were reduced to eps = 0.0025 m (threshold water depth above which cells are wet), hmin = 0.1 m (threshold water depth for Stokes drift), hswitch = 0.05 m (switched from wet to dry) and dzmax = 0.05 m/s/m (maximum bed-level change due to avalanching).
The morphodynamic calibration was performed by comparing the final dune profiles and the percentages of deviation in post-dune volume per meter (ΔV), which was The comparison with the upscaled measured H rms,hf wave heights shows good overall agreement for each hydrodynamic condition, with the bias varying between −0.3 cm (DC) and 5.77 cm (SO). The most agreement with the observation was achieved for the overwash cases with R 2 = 0.78 for SO and R 2 = 0.82 for DO. With respect to nearshore wave heights, the model underestimated the maximum observed wave height, except for SO. Overall, the results are very well in line with the model's small-scale application [32], as similar statistics were derived (see Table 2), which proves that the translation to a larger scale is valid with respect to hydrodynamics. Table 2. Comparison of wave statistics between the small-scale (adapted from [32]) and largescale (this study) applications based on the coefficient of determination (R 2 ) and the bias for each hydrodynamic condition.

Morphodynamic Calibration
The morphodynamic calibration was performed based on the upscaled control cases and by reusing the boundary conditions generated during the hydrodynamic calibration (wbctype = reuse). Due to the similar model performance in small-scale and large-scale situations in terms of hydrodynamics, the final morphodynamic parameter settings of Schweiger and Schuettrumpf [32] were used as initial settings in this study; and the critical avalanching slopes were set to wetslp = 0.2 and dryslp = 0.7. However, the depthscale parameter was set to SF = 2 to be consistent with the upscaled model setup (1:2). As a result, the effected morphodynamic parameters were reduced to eps = 0.0025 m (threshold water depth above which cells are wet), hmin = 0.1 m (threshold water depth for Stokes drift), hswitch = 0.05 m (switched from wet to dry) and dzmax = 0.05 m/s/m (maximum bed-level change due to avalanching).
The morphodynamic calibration was performed by comparing the final dune profiles and the percentages of deviation in post-dune volume per meter (∆V), which was calculated for both measurements and the simulation for the area above the concrete wall (zb > 0.54 m). Furthermore, the Brier skill score (BSS) was calculated: with zb c being the computed bed level, zb m the measured bed level and zb 0 the initial bed level. The BSS defines a model skill as bad (<0), poor (0-0.3), fair (0.3-0.6), good (0.6-0.8) or excellent (0.8-1) according to the classification of van Rijn et al. [48]. Figure 4 shows the initial (dotted) and the final measured (dashed) dune profiles for each hydrodynamic condition (a-d) in combination with various simulation results (colored dash-dotted). Applying the initial settings (blue lines) led to an overprediction of the observed dune erosion for each hydrodynamic condition which was greater for the conditions targeting collision (Figure 4a,c). Similarly to the small-scale application, however, no overwash was computed for each hydrodynamic condition, which is contrary to initial expectations due to the large-scale application. Since the aim of this study was to validate the root model for large-scale and overwash in particular, we chose to induce the occurrence of overwash by additionally varying hydrodynamic parameters, although the model was already hydrodynamically calibrated. calculated for both measurements and the simulation for the area above the concrete wall (zb > 0.54 m). Furthermore, the Brier skill score (BSS) was calculated: with zbc being the computed bed level, zbm the measured bed level and zb0 the initial bed level. The BSS defines a model skill as bad (<0), poor (0-0.3), fair (0.3-0.6), good (0.6-0.8) or excellent (0.8-1) according to the classification of van Rijn et al. [48]. Figure 4 shows the initial (dotted) and the final measured (dashed) dune profiles for each hydrodynamic condition (a-d) in combination with various simulation results (colored dash-dotted). Applying the initial settings (blue lines) led to an overprediction of the observed dune erosion for each hydrodynamic condition which was greater for the conditions targeting collision (Figure 4a,c). Similarly to the small-scale application, however, no overwash was computed for each hydrodynamic condition, which is contrary to initial expectations due to the large-scale application. Since the aim of this study was to validate the root model for large-scale and overwash in particular, we chose to induce the occurrence of overwash by additionally varying hydrodynamic parameters, although the model was already hydrodynamically calibrated. Therefore, various hydrodynamic parameters were additionally considered in order to induce overwash. The majority of these parameters (e.g., alpha (wave dissipation coefficient): 0.5-2; gammax (maximum ratio wave height to water depth): 2, 3; n (power in Roelvink dissipation model): 10,15), and other processes (swrunup (short-wave runup): 1 with facrun = 2) had a negligible effect on the results. However, increasing the delta (δ) parameter resulted in the occurrence of overwash for SO and DO. This parameter, which is zero by default, defines the fraction of Hrms wave height that is added to the water depth h within the calculation of the maximum wave height Hmax: where γ is the breaker index. While overwash was observed for delta ≥ 0.5, a final value of delta = 0.7 was chosen due to the highest agreement with the post-observed dune crest Therefore, various hydrodynamic parameters were additionally considered in order to induce overwash. The majority of these parameters (e.g., alpha (wave dissipation coefficient): 0.5-2; gammax (maximum ratio wave height to water depth): 2, 3; n (power in Roelvink dissipation model): 10,15), and other processes (swrunup (short-wave runup): 1 with facrun = 2) had a negligible effect on the results. However, increasing the delta (δ) parameter resulted in the occurrence of overwash for SO and DO. This parameter, which is zero by default, defines the fraction of H rms wave height that is added to the water depth h within the calculation of the maximum wave height H max : where γ is the breaker index. While overwash was observed for delta ≥ 0.5, a final value of delta = 0.7 was chosen due to the highest agreement with the post-observed dune crest height for SO. In addition, applying delta = 0.7 reduced the erosion at the dune front for all hydrodynamic conditions (Figure 4a-d, yellow). However, the downside of adjusting the delta parameter is a decrease in model accuracy concerning nearshore wave heights, which is reflected by the wave statistics in Table 3. In favor of accurately modeling the relevant overwash conditions, we chose to use a value of delta = 0.7 in the remainder of the investigation. The reason for this was that the validation of the root model for a largescale application was based on the comparison of the belowground with the respective control cases. Thus, deviations in terms of hydrodynamics were assumed to apply to both cases equally. In order to further decrease the erosion at the dune front, the facua parameter was increased to 0.3 (Figure 4, orange lines), which promotes the effect of wave skewness and asymmetry on the sediment advection velocity, and hence increases the onshore-directed sediment transport [49]. However, there was still too much erosion at the dune front for each hydrodynamic condition. Therefore, the increase to facua = 0.3 was combined with delta = 0.7, which further improved the results (Figure 4, purple). For the C-conditions, erosion at the dune front decreased and the computed dune face was in line with the observation for SC1. With respect to statistics (see Table 4), however, model skill remained low (BSS SC1 < −1, BSS DC = −0.24), although the deviation in post-volume decreased considerably from −152% to −32.9% (SC1) and from −98% to 4.1% (DC). For the conditions targeting overwash, the computed post-dune crest heights are in good agreement with the measurements. Nevertheless, low skill was achieved for SO (BSS < −1), although the deviation in post-volume decreased from −78.9% to −4.5%. For DO, however, the final dune profile shape is very well in line with the observation, which demonstrates excellent model skill (BSS = 0.88). The deviation in post-volume (14.6%) can be attributed to the general model-data mismatch between 3 m < x < 9 m (d).

Modeling Aboveground and Belowground Vegetation Cases
In order to model the effect of BGB on dune erosion volumes at large scales, the upscaled belowground cases of Bryant et al. [34] served as comparative examples, where uniformly integrated coconut husk fibers represented belowground biomass. Within the scope of this study, it was hence assumed that the influence of these fibers on the erosion resistance of the soil scales linearly with the model setup. Therefore, it was assumed that the comparison of the simulation results with the upscaled measured final dune profiles (BGB) would be valid, as would the comparison of the dune profile differences (BGB to control).
Due to the inconsistent model performance during calibration in the collision and overwash conditions, the modelling of the belowground cases was performed separately. For the C-conditions, only the constant (with respect to rcc value) root model was applied with a constant and a locally varying rooting depth (z root ). Due to the non-occurrence of overwash in the control simulations, locally increased friction coefficients were not considered. For the O-conditions, however, both the constant root model (with constant/varying z root ) and the friction approach were considered in the simulation. In general, the root model was applied with a maximum rooting depth of z root = 1.15 m, which equals the difference between the dune crest and the height of the sloping wall. In the case of local variation, the rooting depth was defined as the difference between the initial dune profile and the height of the concrete wall (zb = 0.54 m), which was provided for each grid cell through an external file (see Figure 1).
Since locally increasing the bed friction coefficients is the common practice for the consideration of aboveground vegetation in (XBeach) modeling [22,[27][28][29]50], the aboveground cases (AGB) of Bryant et al. [34] served as an additional source of comparison when solely using the friction approach. In the physical experiment, aboveground biomass was represented by 30.5 cm long wooden dowels which were buried half of their length deep. Finally, the combination of the root model and locally increased friction coefficients was compared to the upscaled above and belowground cases (ABGB) of Bryant et al. [34]. In this regard, it should be noted that there was no physical connection between the wooden dowels (ABG) and the husk fibers (BGB) in the physical experiment.

Collision Regime
Applying the root model to the belowground cases of the conditions targeting collision decreased the dune erosion at the dune front for both SC1 and DC. However, the influence of the root model remained constant for rcc ≥ 0.75. As an example, Figure 5 shows the final dune profiles of the measurements (black), the control simulation (blue) and when applying the root model with rcc = 0.75 and a constant (orange) or a spatially varying rooting depth (yellow) for SC1 (a) and DC (b). Given as well are the dune profile differences between belowground and control for SC1 (c) and DC (d). In the case of SC1, using the root model with constant z root resulted in a similar maximum decrease of ∆zb max = 30 cm as the measurement. However, the location of maximum decrease was shifted offshore by one meter and erosion at the dune front (8 m < x < 10 m) was too high in comparison to the upscaled measurement (c). The use of spatially varying rooting depths to decrease the influence of the root model at the dune face had only a minor effect with respect to the control. For DC, similar results can be observed. However, the model predicted a maximum decrease of ∆zb max = 20 cm, which is much lower than the upscaled measurement with ∆zb max = 40 cm (d).

Overwash Regime
Due to the occurrence of overwash in the control simulations, the use of locally increased friction coefficients was considered in the modeling of the O-conditions. In this regard, the friction coefficients were locally increased between 0.035 s/m 1/3 and 0.08 s/m 1/3 (shrubs [28]). Figure 6 shows the simulation results, including the final dune profiles for SO (a) and DO (b) and the dune profile differences (c, d). In general, increasing the bed friction coefficients decreases the erosion on the lee side of the dune. While the measured maximum decrease in erosion (x = 11.4 m) was underestimated by 15 cm for SO when using nveg = 0.08 s/m 1/3 (c, green), the maximum decrease of Δzbmax = 37 cm for DO (d, green) is in good agreement with the upscaled measurement (Δzbmax = 41 cm), but shifted offshore by 1.25 m.
With respect to AGB, the maximum decrease in erosion at x = 11.4 m was underestimated by approximately 5 cm when applying nveg ≥ 0.06 s/m 1/3 (green and purple). For DO, however, the effect of aboveground biomass on dune erosion volumes was overestimated for nveg ≥ 0.035 s/m 1/3 .

Overwash Regime
Due to the occurrence of overwash in the control simulations, the use of locally increased friction coefficients was considered in the modeling of the O-conditions. In this regard, the friction coefficients were locally increased between 0.035 s/m 1/3 and 0.08 s/m 1/3 (shrubs [28]). Figure 6 shows the simulation results, including the final dune profiles for SO (a) and DO (b) and the dune profile differences (c, d). In general, increasing the bed friction coefficients decreases the erosion on the lee side of the dune. While the measured maximum decrease in erosion (x = 11.4 m) was underestimated by 15 cm for SO when using n veg = 0.08 s/m 1/3 (c, green), the maximum decrease of ∆zb max = 37 cm for DO (d, green) is in good agreement with the upscaled measurement (∆zb max = 41 cm), but shifted offshore by 1.25 m.
With respect to AGB, the maximum decrease in erosion at x = 11.4 m was underestimated by approximately 5 cm when applying n veg ≥ 0.06 s/m 1/3 (green and purple). For DO, however, the effect of aboveground biomass on dune erosion volumes was overestimated for n veg ≥ 0.035 s/m 1/3 .

Figure 6. Modeling overwash conditions with belowground biomass: Comparison of initial (thin dotted), final measured control (thick dotted, black) and final measured BG (dashed, black) and AG (dash-dotted, black) dune profiles with simulation results (dash-dotted), including the control simulation (blue) and the application of spatially increased friction coefficients (color) for shallow overwash (a) and deep overwash (b). Profile differences (belowground to control) are shown for the measurements (dashed) and the simulations (dash-dotted) for shallow (c) and deep collisions (d). Note:
Overlay of the purple and green lines. Applying the root model led to results that differed with the overwash conditions and with the use of constant or spatially varying rooting depth. Figure 7 summarizes different final dune profiles (a, b) and dune profile differences (c, d) for different root model configurations compared to the upscaled measured data. For SO, using rcc values of 0.75 and 1 resulted in good agreement between the model and the measurements with respect to the maximum decrease (x = 11.5 m) in the presence of BGB (c). Furthermore, incorporating a spatially varying rooting depth decreased the effect of the root model on the model predicted erosion at the dune front (x < 11 m) in comparison to the simulations with a constant rooting depth (orange, purple). As a result, greater agreement with the measurements was achieved. For the back of the dune (x > 12 m), however, the predicted decrease in dune erosion was lower than the measurement. For DO, the best agreement with respect to dune profile differences between BGB and control was achieved for rcc = 0. 35 (d). Although the location of maximum decrease was shifted by approximately 1.5 m offshore, the general profile shape was in good agreement with the observations of the front and the back of the dune. Contrary to the application to SO, the use of a spatially varying rooting depth led to only minor differences in the dune foot area (8 m < x < 9 m). Applying the root model led to results that differed with the overwash conditions and with the use of constant or spatially varying rooting depth. Figure 7 summarizes different final dune profiles (a, b) and dune profile differences (c, d) for different root model configurations compared to the upscaled measured data. For SO, using rcc values of 0.75 and 1 resulted in good agreement between the model and the measurements with respect to the maximum decrease (x = 11.5 m) in the presence of BGB (c). Furthermore, incorporating a spatially varying rooting depth decreased the effect of the root model on the model predicted erosion at the dune front (x < 11 m) in comparison to the simulations with a constant rooting depth (orange, purple). As a result, greater agreement with the measurements was achieved. For the back of the dune (x > 12 m), however, the predicted decrease in dune erosion was lower than the measurement. For DO, the best agreement with respect to dune profile differences between BGB and control was achieved for rcc = 0. 35 (d).
Although the location of maximum decrease was shifted by approximately 1.5 m offshore, the general profile shape was in good agreement with the observations of the front and the back of the dune. Contrary to the application to SO, the use of a spatially varying rooting depth led to only minor differences in the dune foot area (8 m < x < 9 m). Despite the comparison of dune profile differences, the overwash sediment volumes per meter width (Vover) were calculated for each simulation: Subsequently, the percentage decrease in overwash volume in each control simulation was derived to enable a comparison to the results of the small-scale physical experiment [34]. It should be noted that the latter were obtained by collecting and weighing the overwash sediment and are hence based on total sediment mass. The ratios between BGB and control are summarized in Table 5 for shallow-overwash and in Table 6 for deep-overwash. For shallow-overwash, the presence of ABG and BGB reduced the overwash volume by 34% (37.6 kg to 24.9 kg) and by 58% (37.6 kg to 15.9 kg) in the physical experiment, respectively. The best agreement with the BGB measurement was achieved when applying the root model with rcc = 0.75 with a percentage decrease of 59% (const. zroot) and 60% (var. zroot). In case of the former, this corresponds to a decrease from 0.58 m 3 /m to 0.24 m 3 /m. Hypothetically, the observed decrease in overwash sediment mass of (37.6-15.9) kg = 21.7 kg would be equivalent to an upscaled decrease in overwash volume of 0.04 m 3 /m (the translated overwash volume follows from Vover = 21.7 kg/2650 kg/m 3 /1.5 m⸱7.5 = 0.041 m 3 /m) when assuming a sediment density of 2650 kg/m 3 (XBeach default), and when using the small-scale flume width of 1.5 m and the scale factor of 7.5. Thus, the model theoretically overpredicts the decrease in overwash volume by a factor of 8.5 (= (0.58−0.24)/0.04). For AGB, the use of locally increased friction coefficients underestimates the measured percentage decrease in overwash volume by at least 5% (nveg ≥ 0.06 s/m 1/3 ).
For DO, the presence of AGB decreased the total overwash sediment mass by 20% from 144.2 kg to 115.5 kg. For BGB, the total overwash sediment decreased to 78 kg (46%). As an example, the theoretical upscaling for BGB resulted in an overwash volume decrease of 0.12 m 3 /m for the measurement. With respect to modeling, the greatest Despite the comparison of dune profile differences, the overwash sediment volumes per meter width (V over ) were calculated for each simulation: Subsequently, the percentage decrease in overwash volume in each control simulation was derived to enable a comparison to the results of the small-scale physical experiment [34]. It should be noted that the latter were obtained by collecting and weighing the overwash sediment and are hence based on total sediment mass.
The ratios between BGB and control are summarized in Table 5 for shallow-overwash and in Table 6 for deep-overwash. For shallow-overwash, the presence of ABG and BGB reduced the overwash volume by 34% (37.6 kg to 24.9 kg) and by 58% (37.6 kg to 15.9 kg) in the physical experiment, respectively. The best agreement with the BGB measurement was achieved when applying the root model with rcc = 0.75 with a percentage decrease of 59% (const. z root ) and 60% (var. z root ). In case of the former, this corresponds to a decrease from 0.58 m 3 /m to 0.24 m 3 /m. Hypothetically, the observed decrease in overwash sediment mass of (37.6-15.9) kg = 21.7 kg would be equivalent to an upscaled decrease in overwash volume of 0.04 m 3 /m (the translated overwash volume follows from V over = 21.7 kg/2650 kg/m 3 /1.5 m · 7.5 = 0.041 m 3 /m) when assuming a sediment density of 2650 kg/m 3 (XBeach default), and when using the small-scale flume width of 1.5 m and the scale factor of 7.5. Thus, the model theoretically overpredicts the decrease in overwash volume by a factor of 8.5 (= (0.58−0.24)/0.04). For AGB, the use of locally increased friction coefficients underestimates the measured percentage decrease in overwash volume by at least 5% (n veg ≥ 0.06 s/m 1/3 ). * Please note that these values equal the percentage differences in sediment loss (kg) between belowground and control values within the scope of the original (small-scale) physical experiment (adapted from [34]). For DO, the presence of AGB decreased the total overwash sediment mass by 20% from 144.2 kg to 115.5 kg. For BGB, the total overwash sediment decreased to 78 kg (46%). As an example, the theoretical upscaling for BGB resulted in an overwash volume decrease of 0.12 m 3 /m for the measurement. With respect to modeling, the greatest agreement was achieved when applying the root model with rcc = 0.35 and a variable rooting depth. Then, the predicted overwash volume decreased by 43% from 2.49 m 3 /m to 1.41 m 3 /m. The theoretical comparison between the model and the upscaled measurements found overestimation by a factor of 9 (= (2.49−1.41)/0.12). Concerning locally increased friction, using values of n ≥ 0.06 s/m 1/3 led to a decrease in overwash volume by 41%, and hence, similar (good) agreement with the measurements.
Further analysis was conducted by evaluating the root model (rcc = 0.35, var. z root ) and using locally increased friction coefficients (n veg = 0.06 s/m 1/3 ) for SO after each wave burst. Figure 8 presents the dune profiles (a−c), the cumulative erosion and the bed-level changes due to avalanching (d−f) and the maximum sediment concentration during each wave burst (g−i). After the first wave burst, erosion at the dune crest (x = 11.5 m) was greatest for the control simulation (blue) with a decrease of ∆zb wb1,ctrl = −8 cm (d). Locally increasing the bed friction (purple) reduced erosion at the dune crest to ∆zb wb1,fric = −5 cm. With the root model (green), erosion at the dune crest was almost prevented (∆zb wb1,rm < −1 cm). As a result, the cumulative erosion on the lee side of the dune was greatest for the control simulation with a maximum bed level change of ∆zb wb1,ctrl = −12 cm at x = 12 m, which decreased to ∆zb wb1,fric = −10 cm (friction) or ∆zb wb1,rm = −7 cm (root model). The results show further that there is no contribution of avalanching for x > 11.2 m (dzav = 0). With respect to the maximum sediment concentration, only small differences can be observed between the control and the friction approach. Applying the root model reduced the maximum sediment concentration at both the front (9 m < x < 11.5 m) and the back of the dune.
In the remainder of the control simulation, erosion at the dune crest and on the lee side of the dune increased to a maximum of ∆zb wb2,ctrl = −21 cm after the second wave burst (e) and ∆zb wb3,ctrl = −28 cm after the third wave burst (f). For locally increased friction and the root model, the maximum erosion decreased to ∆zb wb2,fric = −15 cm and ∆zb wb3,fric = −20 cm (friction), and to ∆zb wb2,rm = −10 cm and ∆zb wb3,rm = −12 cm (root model). In terms of maximum sediment concentrations, similar results as for the first wave burst were obtained when using locally increased friction (similar to control) and the root model (lower than control), which also applies to the non-occurrence of avalanching for x > 11.2 m. In summary, using the root model and thus increasing the critical velocity for erosion along the dune reduced the overall erosion, especially at the dune crest. As a result, erosion on the lee side of the dune was reduced during the second and third wave bursts without bed level changes due to avalanching, which finally resulted in lower overwash sediment volumes (see Table 5).
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 15 of 21 bed level changes due to avalanching, which finally resulted in lower overwash sediment volumes (see Table 5).

Modeling the Combination of Aboveground and Belowground Vegetation
Due to the occurrence of overwash and thus an influence of locally increased friction coefficients on the predicted erosion volumes, further simulations were conducted by using both the root model and locally increased friction coefficients. In this way, the root model was applied with varying rooting depth and with rcc = 0.75 for SO and rcc = 0.35 for DO based on the BGB modeling (see Figure 7). The results were compared to the upscaled ABGB results of Bryant et al. [34]. Figure 9 presents the final dune profiles for SO (a) and DO (b) and the dune profile differences between ABGB and the control (dash-dotted), and between the simulation and the measurements (dotted in (c) and (d)). Furthermore, the percentage decreases in sediment loss to overwash are summarized in Table 7. In general, the comparison of final dune profiles (a, b) shows that combining the two approaches decreased dune erosion volumes on the lee side of the dune compared to using the root model alone (orange). For SO, the combination of the root model with locally increased friction coefficients to nveg = 0.08 s/m 1/3 underpredicted the upscaled measured decrease due to the presence of ABGB (c). Nevertheless, the final dune profile is in good agreement with the measurement for x > 11.3 m, as shown by the differences between simulation and measurement (dotted line in c). However, the observed decrease in overwash volume (86%) was underpredicted by 11%. Figure 8. Time evolution of various output parameters for SO including dune profiles after the first (a), second (b) and third (c) wave burst, the cumulative erosion and the bed level change due to avalanching after each wave burst (d-f) and the maximum sediment concentration during each respective wave burst (g-i).

Modeling the Combination of Aboveground and Belowground Vegetation
Due to the occurrence of overwash and thus an influence of locally increased friction coefficients on the predicted erosion volumes, further simulations were conducted by using both the root model and locally increased friction coefficients. In this way, the root model was applied with varying rooting depth and with rcc = 0.75 for SO and rcc = 0.35 for DO based on the BGB modeling (see Figure 7). The results were compared to the upscaled ABGB results of Bryant et al. [34]. Figure 9 presents the final dune profiles for SO (a) and DO (b) and the dune profile differences between ABGB and the control (dash-dotted), and between the simulation and the measurements (dotted in (c) and (d)). Furthermore, the percentage decreases in sediment loss to overwash are summarized in Table 7. In general, the comparison of final dune profiles (a, b) shows that combining the two approaches decreased dune erosion volumes on the lee side of the dune compared to using the root model alone (orange). For SO, the combination of the root model with locally increased friction coefficients to n veg = 0.08 s/m 1/3 underpredicted the upscaled measured decrease due to the presence of ABGB (c). Nevertheless, the final dune profile is in good agreement with the measurement for x > 11.3 m, as shown by the differences between simulation and measurement (dotted line in c). However, the observed decrease in overwash volume (86%) was underpredicted by 11%. thermore, only minor errors occurred in the simulations using both the root model and the friction approach. This also applies to the percentage decrease in sediment loss to overwash, which was 20% (nveg = 0.035 s/m 1/3 ) to 18% (nveg = 0.08 s/m 1/3 ) lower than observed in the physical experiment (83%). In order to increase the effect of the root model, combining rcc = 0.5 with a local friction value of nveg = 0.08 s/m 1/3 resulted in good agreement with respect to dune profile differences (d) and increased the percentage decrease in overwash volume to 70.4%. Figure 9. Modeling overwash conditions with aboveground and belowground biomass: Comparison of initial (thin dotted), final measured control (thick dotted, black) and final measured ABGB (dashed, black) dune profiles with simulation results (dash-dotted, black) for shallow-overwash (a) and deep-overwash (b). The root model was applied with a varying rooting depth (varzroot) and combined with locally increased friction coefficients. Further shown are the profile differences between belowground and the control for the measurements (dashed) and the simulations (dash-dotted), and the difference between simulations and measurements (dotted) for shallow-(c) and deep-collisions (d). Note: Overlay of the green and purple lines.

Discussion and Conclusions
Due to a general model-data mismatch in small-scale applications [32], this study aimed to validate the root model at larger scale based on an upscaled model setup of Figure 9. Modeling overwash conditions with aboveground and belowground biomass: Comparison of initial (thin dotted), final measured control (thick dotted, black) and final measured ABGB (dashed, black) dune profiles with simulation results (dash-dotted, black) for shallow-overwash (a) and deep-overwash (b). The root model was applied with a varying rooting depth (varzroot) and combined with locally increased friction coefficients. Further shown are the profile differences between belowground and the control for the measurements (dashed) and the simulations (dash-dotted), and the difference between simulations and measurements (dotted) for shallow (c) and deep-collisions (d). Note: Overlay of the green and purple lines. While for DO, the root model resulted in good agreement with the upscaled BGB measurement (see Figure 7), the local increase of friction coefficients overpredicted the effect of AGB with n veg = 0.035 s/m 1/3 (see Figure 6). Concerning ABGB, however, combining the root model (rcc = 0.35) with locally increased friction coefficients underpredicted the measured decrease in dune erosion volume, even for n veg = 0.08 s/m 1/3 (Figure 9d). Furthermore, only minor errors occurred in the simulations using both the root model and the friction approach. This also applies to the percentage decrease in sediment loss to overwash, which was 20% (n veg = 0.035 s/m 1/3 ) to 18% (n veg = 0.08 s/m 1/3 ) lower than observed in the physical experiment (83%). In order to increase the effect of the root model, combining rcc = 0.5 with a local friction value of n veg = 0.08 s/m 1/3 resulted in good agreement with respect to dune profile differences (d) and increased the percentage decrease in overwash volume to 70.4%.

Discussion and Conclusions
Due to a general model-data mismatch in small-scale applications [32], this study aimed to validate the root model at larger scale based on an upscaled model setup of Bryant et al. [34]. It was expected that the overall performance and especially the predictions of overwash would improve, since XBeach's default parameters are optimized for large-scale modeling. With respect to H rms,hf wave heights, good agreement was achieved between the model and (upscaled) measurements for each hydrodynamic condition. Similarly to the small-scale application, dune overwash was not computed regardless of the hydrodynamic condition, which did not change when activating short-wave runup. Only the increase of the delta parameter, and thus, the maximum wave height in the nearshore zone, resulted in the occurrence of dune overwash for the O-conditions. From this follows that the general model-data mismatch described in Schweiger and Schuettrumpf [32] might not only result from the small-scale application, but might (also) be due to other model inaccuracies, such as an underestimation of wave runup, as discussed in previous studies [51][52][53].
Although adjustments to various model parameters were considered in the morphodynamic calibration, erosion volumes at the font of the dune were significantly overestimated for all conditions except DO (see Figure 4). As a result, low model skill was achieved for SC1, DC and SO (BSS < 0), whereas an excellent skill was achieved for DO (BSS = 0.88). However, the overestimated erosion at the dune face and toe is in line with the findings of van Dongeren et al. [54], who, among others, observed an overprediction of erosion around the mean waterline, possibly due to inaccuracies in the modeling of sediment motion in the swash zone (see also [53,55]). In this study, the overestimation of erosion at the dune toe could have been further reduced by further increasing the facua parameter, which promotes the onshore-directed sediment transport. However, this would have been followed by an increasing deviation of the dune front position between the model and measurements, which was already observed for SO when applying facua = 0.3 with delta = 0.7 (see Figure 4b). Therefore, deviations at the dune front were accepted since the main focus in this study was on the erosion at the dune crest and on the lee side of the dune due to overwash.
The effects of locally increased friction coefficients and the root model on dune erosion volumes at a large scale were validated by using the upscaled biomass cases (AGB, BGB, ABGB) of Bryant et al. [34] as a reference. Therefore, a major assumption in this study was that the effects of the wooden dowels (AGB), the husk fibers (BGB) and their combination (ABGB) on the erodibility of the sandy soil scale linearly with the translation to a large scale. From this follows that the validation of two approaches to consider land-based biomass in this study was not based on real data from a physical experiment but relied on a theoretical comparison. Nevertheless, applying the root model to the O-conditions resulted in good agreement with the upscaled BGB measurements, which was better when spatially varying rooting depths were applied. However, the best results for SO and DO were achieved with different rcc values. While for SO a value of rcc = 0.75 led to good agreement concerning dune profile differences (BGB to control), the application to DO required a reduction to rcc = 0.35. Although these results are based on a theoretical comparison, the effect of the root model on the model predicted erosion should be consistent and independent of the hydrodynamic conditions. Therefore, further investigations are necessary to improve the performance of the root model. Due to a lack of appropriate data concerning the interaction between belowground biomass and sandy soil, future studies should focus on quantifying the effects of land-based (belowground) biomass on the erodibility of a dune system, as suggested by Figlus et al. [33]. With respect to further validating the root model in XBeach, these studies should be conducted at a large scale and additionally consider real dune vegetation.
Using locally increased friction coefficients decreased the dune erosion volumes on the lee side of the dune. Concerning BGB modeling, using a local friction coefficient of n veg = 0.08 s/m 1/3 underestimated the decrease in erosion, which was greater for SO (see Figure 6). Since the friction approach is the common practice for considering the effect of land-based vegetation on overwash-dominated erosion in (XBeach) modeling [22,[27][28][29]50], the results were also compared to the AGB cases. For DO, the predicted decrease in erosion was already overestimated for n veg = 0.035 s/m 1/3 . For SO, however, the influence of AGB on dune erosion volumes was underestimated even when using the maximum value (n veg = 0.08 s/m 1/3 ) that was considered in this study. Therefore, the dynamic roughness module [30], which varies the Manning roughness in vegetated areas due to burying and erosion, was not considered in this study. Since the BGB measurements were underestimated when using the friction approach (see Figure 6c,d), the greater effect of ABGB on the dune erosion volumes would presumably have led to an even greater underestimation.
Nevertheless, the combination of the root model with locally increased friction coefficients was validated with the upscaled ABGB measurements. Due to the good performance in the BGB modeling (root model with rcc = 0.35) and the overestimation observed in the AGB modeling for DO (friction approach), combining the two approaches was expected to result in good agreement with the ABGB measurements for DO. Although less erosion was predicted on the lee side of the dune in comparison to the sole use of the root model, the overall decrease due to the presence of ABGB was lower in comparison to the upscaled measurement and less sensitive to variations in local bed friction (see Figure 9d). In this regard, increasing the rcc value to 0.5 resulted in a good agreement between model and upscaled measurements. This necessary increase could be attributed to the general model setup, where the wooden dowels were buried half way. Thus, the belowground parts of the dowels could have further increased the resistance of the dune in addition to the effect of the uniformly integrated coir fibers, hence requiring a higher contribution of the root model, and hence, higher rcc values.
In summary, the application of the XBeach model to an upscaled model setup of Bryant et al. [34] showed a general model-data mismatch for all hydrodynamic conditions except for DO, which was mainly due to overestimated dune erosion volumes around the initial waterline. Nevertheless, applying the root model to the upscaled BGB cases reduced the predicted dune erosion for all hydrodynamic conditions. In this regard, incorporating spatially varying rooting depths further improved the results at the dune front, although the overall effect of the root model differed between the simulations. However, there is a lack of appropriate data to further validate and improve the root model. Although various wave flume experiments have been conducted in recent years, there is a specific need for research focusing on the effect of land-based biomass on the erodibility of dunes during collision and overwash at large scales. In particular, the separate investigation of aboveground and belowground biomass concerning wave-induced dune erosion at large scales and the individual contributions of different plant characteristics (e.g., rooting depth, plant/root density and maturity) would allow one to evaluate the specific influences of these plant-related parameters on the erosive processes. Furthermore, these studies should additionally address the effects of land-based vegetation on the temporal evolution of dune scarping, avalanching and failure.