Salt Marsh Establishment and Eco‐Engineering Effects in Dynamic Estuaries Determined by Species Growth and Mortality

Growth conditions and eco‐engineering effects of vegetation on local conditions in coastal environments have been extensively studied. However, interactions between salt marsh settling, growth, and mortality as a function of hydromorphology and eco‐engineering lack sufficient understanding to forecast morphological development of dynamic systems. We predict salt marsh establishment with an ecomorphodynamic model that accounts for literature‐based seasonal settling and life‐stage‐dependent growth and mortality of a generic salt marsh species. The model was coupled to a calibrated hydromorphodynamic model of an intertidal bar and, on a coarser grid, to the entire Western Scheldt estuary. To quantify the importance of eco‐engineering effects we compared the dynamic model results to a static model approach. The ecomorphodynamic model reproduces spatial pattern, cover, and growth trends over 15 years. The modeled vegetation cover emerges from the combination of a positive and a new negative eco‐engineering effect: vegetation reduces tidal flow strength facilitating plant survival while the developing salt marsh increases the hydroperiod, which limits large‐scale marsh expansion. The reproduced spatial gradient in vegetation density by our model is strongly correlated to their life‐stages, which underlines the importance of age‐dependence when modeling vegetation and for predictions of the stability of the marsh. Upscaling of the model to the entire estuary on a coarser grid gives implications for grid size‐dependent modeling of hydrodynamics and vegetation. In comparison with static model results, the eco‐engineering effects reduce vegetation cover, showing the importance of vegetation dynamics for predictions of salt marsh growth.

Salt marshes maximize their chance of survival by adaptation to their environment across different spatial and temporal scales (Corenblit et al., 2015;Holling, 1973). At the scale of patches, physical plant properties, such as stem rigidity, height, diameter, and density mainly determine the effects of vegetation on flow and sediment transport through flow resistance. This, in turn, alters the environment for the species, which is called eco-engineering effects (Jones et al., 1994). As flow reduction within patches increases the settlement of suspended particles between plant stems and leaves, an increase in bed elevation through sedimentation of suspended organic and mineral sediments is promoted (Reed, 1990). Since increased elevation reduces inundation stress, these processes were previously referred to as local, or small-scale, positive eco-engineering effects (Bouma et al., 2009). At the same time, salt marsh plants can facilitate channel 10.1029/2019JF005092 formation through erosion between salt marsh patches shaping the drainage network of the entire marsh (Schwarz et al., 2015;Temmerman et al., 2007) and can lead to large-scale self-organization of the system ( Van de Koppel et al., 2012). The necessity to model such systems on the reach scale is also demonstrated for fluvial systems that show similar positive and negative eco-engineering effects, wherein vegetated bars and floodplains affect channels and floodplains over much larger distances than the scale of patches through backwater effects Van Oorschot et al., 2016). The intensity of these processes is strongly dependent on the density of the vegetation (Leonard & Luther, 1995;Van Wesenbeeck et al., 2008) and vegetation growth (Bouma et al., 2013). Physical plant properties of salt marshes, located in temperate regions, show seasonal variations in biomass, characterized by seedling germination in spring, peak biomass at the end of summer, and senescence in winter (Ibáñez et al., 2012;Morris et al., 2002). With such dynamic vegetation development, the magnitude of eco-engineering effects will not only change through space but also undergo considerable change throughout the growth season and between years.
Before the eco-engineering effects can occur, initial settling and establishment of salt marsh take place. In particular, seedling survival requires periods with lower disturbance to allow sufficient root growth to withstand pressures such as currents or waves (Balke et al., 2014;Cao et al., 2018). In general, settling, growth, and mortality are species-specific functions of environmental conditions and life-stage-dependent stress tolerance to inundation, desiccation, flow velocity, scour, and burial (Bouma et al., 2013;Friedrichs & Perry, 2001;Fagherazzi & Sun, 2004;Schwarz et al., 2018;Van Hulzen et al., 2007). Seedling survival was shown to be mainly influenced by the combination of burial, erosion, and flow velocity (Wang & Temmerman, 2013;Willemsen et al., 2018), while survival of mature plants was shown to be mainly influenced by inundation time expressed in either flooding frequency or hydroperiod (Balke et al., 2016;D'Alpaos et al., 2006;Hughes et al., 2012;Morris & Haskin, 1990;Morris et al., 2002;Mendelssohn & Morris, 2002;Reed, 1990). Inundation time and mature salt marsh development depend on environmental factors such as daily to fortnightly variations in the tide and annual water level variations (Chapman, 1964;Morris et al., 2002;Mudd et al., 2004;Silvestri & Marani, 2004;Suchrow & Jensen, 2010).
Salt marsh response to future changes likely depends on the dynamics of the morphology. In systems with little morphological change on a timescale of decades to centuries, the potential for a salt marsh to keep up with sealevel rise depends on fine sediment supply (Kirwan et al., 2010). However, the feedback between vegetation dynamics and eco-engineering response remains poorly understood. While the abiotic properties alone could be used to predict locations where salt marshes will establish, this ignores the feedback between the salt marsh and abiotic stressors, that is, the eco-engineering effects. Recent research conducted in the Western Scheldt estuary, the Netherlands, showed that the border between the high biomass high marsh and low biomass pioneers zone can be found at a relative inundation period of around 0.45, which was shown to be applicable to Dutch as well as to North American salt marshes (Van Belzen et al., 2017). This empirical threshold is the result not only of the individual plant properties but also of the modifications of the abiotic conditions by the salt marsh. In other words, the survival of plants depends on abiotic spatial and temporal variables (i.e., water levels, sediment supply, salinity) that in turn are affected by vegetation distribution and collective plant characteristics of the life-stages of all species present Cowles, 1911;Wang & Temmerman, 2013). We aim to disentangle the interaction between vegetation establishment, growth, mortality, on the one hand, and of specific eco-engineering effects on salt marsh establishment, on the other hand. In particular, we study a dynamic estuarine environment where salt marsh established only recently and conditions are well described. This will contribute to a better understanding and system-scale predictability of initial colonization of pioneer vegetation on bare tidal flats and for managed realignment projects.
To address this objective, a numerical model is needed that incorporates the abovementioned interactions between flow, morphology, and vegetation over time and space, further referred to as ecomorphodynamic model. Past work addressed only part of the processes and feedbacks; for instance, some models consider vegetation flow interactions but model spatial vegetation establishment as constant vegetation biomass determined by tidal benchmarks, such as inundation period or bed elevation, (e.g., Mudd et al., 2004;D'Alpaos et al., 2006D'Alpaos et al., , 2007Marani et al., 2010;D'Alpaos, 2011), further referred to as static establishment approach. Others include detailed interactions between vegetation flow and inundation time but oversimplified vegetation-induced sedimentation and erosion (e.g., Rodríguez et al., 2017). In long-term modeling, growth and seasonality are usually simplified, where the salt marsh development is updated yearly or quarterly in, for example, Temmerman et al. (2007)  is introduced accounting for reduced biomass in winter (Best et al., 2018;Mudd et al., 2004). Simplified bathymetries, hydrodynamics, and sediment transport computations limit the development of natural vegetation response in all these models, and the adaptation of the plants resilience with aging from seedling to mature life-stages has not been accounted for in any ecomorphodynamic model so far.
Here we apply a new ecomorphodynamic model framework, able to simulate spatiotemporal salt marsh development based on detailed interactions with hydrodynamics, sediment transport, field morphology, and field forcings. Vegetation parameters, for example, growth and mortality, are literature-based, season-dependent and age-or life-stage-dependent, which will allow local and far-field eco-engineering effects to emerge from frequent coupling to the hydromorphodynamics. We test the ability of this new ecomorphodynamic model to predict salt marsh development on both an intertidal bar and on the entire Western Scheldt estuary, for which detailed field data for validation are available.

Methods
To determine the main processes that govern salt marsh establishment in dynamic systems we investigated an intertidal bar, located in the Dutch part of the Western Scheldt estuary, southwest of the Netherlands, as a reference case. For the Western Scheldt, a wide range of monitoring data (bed elevations and vegetation distribution) are available, allowing validation of model results for several years of salt marsh establishment and growth including the feedbacks with the morphological development.
We developed a new ecomorphodynamic model consisting of a dynamic vegetation model, which is coupled with a hydromorphodynamic model (HM-model). The dynamic vegetation model was originally developed by Van Oorschot et al. (2017), as a dynamic model for riparian trees and herbs. Here we extended the dynamic vegetation model for generic salt marsh growth and mortality for two life-stages and implemented a coupling interval with the HM-model every tidal cycle of the dominant M2-tide. Every coupling, the hydromorphodynamic calculations were used as input for the dynamic vegetation model, which updated the spatial vegetation distribution. The new vegetation distribution was fed back into the HM-model as a hydraulic roughness and an extra drag term. As a result, we captured the feedback loop between vegetation dynamics and morphology on dense temporal scales that allowed for a detailed study of the co-occurring processes.

Site Description
The investigated intertidal flat is a 4 km long bar located in the 160 km long Western Scheldt estuary, the Netherlands (Figure 1). The dynamic tidal bar of Walsoorden experienced recent salt marsh establishment In (a) the original morphological classes are displayed with a distinction between high and low energetic environments, littoral classes along bed elevation, and salt marsh cover. Based on these classes we determine classes that included either vegetation or bare soil with bare soil being further classified by sediment type into either muddy sediments or sandy sediments (b). For our analysis we compare our results to (b). starting from the 1990s and has been intensively studied for understanding its bed form patterns, morphological, habitat and salt marsh development, and macrobenthic dynamics (Cleveringa, 2014;De Vet et al., 2017;Plancke et al., 2010;Van der Wal et al., 2011.
The Western Scheldt estuary is a tide-dominated, mesotidal to macrotidal estuary with an upstream increasing tidal range from 3.8 m at the mouth to 5.2 m near the city of Antwerp (Meire et al., 2005). The importance of wave forcing on morphology is gradually exceeded by the effects of tidal currents from the mouth toward the mesohaline part of the estuary (Hu et al., 2018). The estuary has been impacted by dredging and dumping since the 1970s (Van Damme et al., 1999). The dredging activity in combination with local sea level rise is thought to have caused higher bars and steeper bar margins and caused the intertidal zone in the estuary to constantly decrease in dimensions (De Vriend et al., 2011;van Dijk, Cox, et al., 2019). The marshes in the Western Scheldt persisted because of their fast adaptation to changes in the hydrological regime (Temmerman, Bouma, Govers, & Lauwaet, 2005;Van Damme et al., 2001;Wang & Temmerman, 2013

10.1029/2019JF005092
The study site, the tidal bar of Walsoorden, is located in the mesohaline part of the estuary around 50 km land inward from the mouth between an ebb and flood channel of the Western Scheldt estuary (Figure 1). It consists mainly of sandy sediments of 50-150 μm with less than 10% mud content (Van Eck, 1999). Ecotope maps were collected on a multiannual basis to monitor locations of salt marshes (between high water and spring tide), bare tidal flats (between low and high water), and open water (shallow and deep) based on aerial imagery and relative bed elevation (Van Damme et al., 2001). For our analysis and comparison to the model results we simplified our map to the relevant morphological classes present in the ecotope maps (Figure 2), which are bare soil (with one class for sediments classified as muddy), "sparse/pioneer" vegetation and "dense/mature" vegetation distinguished by high and low coverage (dense: >50% and sparse: <50%).
In 1996, patches of the pioneer species Spartina anglica and Salicornia ssp. were observed for the first time by Stikvoort (1996) and over time led to the development of a mature salt marsh on the bar of Walsoorden (Cleveringa, 2014). The strongest increase in pioneer vegetation occurred between 2004 and 2008, which coincides with reduced high water levels between 2001 and 2005 (Balke et al., 2014). At the same time, dredging and dumping were started, nourishing the area around the bar with little effect on ecology (Ides et al., 2007;Van Der Wal et al., 2005, 2011. Overall, a steepening of the bar margins with an increase in average bed elevation since 2001 was observed (Cleveringa, 2014;De Vet et al., 2017). This combined development of bed accretion and increase in vegetation cover gives the opportunity to gain insights into the factors that control successful salt marsh establishment of pioneer species characteristic for NWEurope.

The Ecomorphodynamic Model Setup
Our dynamic vegetation model was coupled to a depth-averaged two-dimensional hydromorphological model in Delft3D (version 4.01.00), a code that solves the shallow water equations (Lesser et al., 2004;van Dijk, Hiatt, et al., 2019). As sediment transport predictor we used Van Rijn (2004) as implemented in Delft3D ( van Dijk, Hiatt, et al., 2019). The morphodynamics are included by solving the advection-diffusion equation for suspended sediment based on reference concentration and capacity transport for bed load combined with mass conservation (see Lesser et al., 2004 for more detail). The interaction between vegetation and flow is realized by the existing trachytope module in Delft3D, which calculates an average flow resistance in each grid cell from the present fractions of vegetation, and their properties, based on the Baptist equation (Baptist et al., 2007).
The subdomain of the hydromorphodynamic (HM) model is taken from a model of the Western Scheldt (Nederlands-Vlaams [NeVla] or Dutch-Flemish model) that was implemented in Delft3D, hydrodynamically calibrated (Maximova, Ides, Vanlede, et al., 2009;Vroom et al., 2015) and later optimized for morphology (Grasmeijer et al., 2013;Schrijvershof & Vroom, 2016). For investigating the intertidal bar at a higher grid resolution, the hydromorphodynamic domain was decomposed into two domains to reduce computational time while allowing for accurate hydrodynamic and vegetation computations. The outer domain comprises of grid cell size up to 180 m. We applied the vegetation modeling on the area of interest on a refined inner grid (down to 16 m cell size) that covers the entire shoal and partly the adjacent channels (see Figure 1). Our hydrodynamic time step was 6 s to allow for stable simulations. To test the possibility for upscaling we also used the original coarser, computationally cheaper, NeVla-model grid for the region from the mouth to the Dutch-Belgian border (see blue-red area in Figure 1).
The boundary conditions of the outer domain were time series of seaward water level measurements and landward discharge measurements in the Western Scheldt in 2013 ( van Dijk, Hiatt, et al., 2019). Consequently, most tidal constituents, wind and wave surges, and seasonal discharge variations were captured in the boundaries. However, waves were not solved in the model as the Western Scheldt is a tide-dominated system (Davis Jr & Hayes, 1984) and a small fetch leads to limited wave effects on the tidal bars (Hu et al., 2018). To solve the problem of fitting the tidal period in an integer number of days (Duran-Matute & Gerkema, 2015) we compressed the period of the dominant M2-tide to 12 hr, which led to a reduced forcing period of each water level by 3%. That way we were able to couple our vegetation model each tidal signal, leading to an ecological time step of 12 hr of hydrodynamics. This step was an important prerequisite to provide continuous forcing for each vegetation computation that included the entire variation along one M2-tide and at the same time allowed to upscale the model to morphological timescales.
From the boundaries we chose four representative spring-neap cycles (i.e., MHW, MW, and MLW are in the same range as for the entire year) that we upscaled to four morphological years using a morphological acceleration factor of 24. The morphological acceleration factor multiplies the morphological computations by a 10.1029/2019JF005092 factor to accelerate morphological development compared to the hydrodynamic forcing (Lesser et al., 2004). As a result, we assume that one tidal signal (here now 12 hr) roughly represents one spring-neap cycle in morphological time (12 days). We apply 28 tidal signals per morphological year to fit the water levels of an entire spring-neap cycle into each year (14 days of hydrodynamic boundaries). Consequently, one morphological year is represented by one spring-neap cycle of hydrodynamic boundaries. As morphological and ecological time were directly coupled, this resulted in 28 vegetation updates each morphological year. The four spring-neap cycles led into four morphological, or ecological, years for each scenario with a total simulation time of 56 days of hydrodynamic boundaries. For our longer calculations we repeated the prescribed boundary conditions and ran 15 morphological years, which was similar to the period for which field data were available.
The hydrodynamic model computes water levels and flow velocities through mass and momentum conservation. With each ecological update, the settlement locations and surface fractions of the vegetation and the physical properties (height, stem diameter, and density) were fed back into the HM-model. All life-stages in a grid cell were combined into a single Chézy coefficient C for each cell and used for the calculation of an additional drag term 2 u 2 (m/s 2 ) in the momentum equation: (1/m) is calculated based on the relative change in velocity, which depends on the water depth h (m), vegetation height h v (m) and net bed roughness C ( √ m/s): where C b ( √ m/s) is real bed roughness, C D (-) the bulk drag coefficient accounting for leaves and branches of the plant, and n (m/m 2 ) plant density multiplied by stem diameter. C is defined for the case of emerged or submerged vegetation as with g = 9.81 (kg/m 2 ) gravity and = 0.41 (-) von-Kármán constant. For simplicity we chose 1 and 1.1 for C D for juveniles and mature life-stage, respectively, approximating the plants as smooth rigid cylinders (Baptist et al., 2007). From equations (2a) and (3a) it is obvious that the seasonal changes of plant height and stem diameter allow for a detailed temporal and spatial representation of the vegetation growth that to-date has not been included in numerical modeling of salt marsh vegetation.
To compute the total vegetation effect in each cell, both C and are weighted by the fraction of the cell that is covered by the vegetation. For the computation of several vegetation fractions within one cell the average weighted values are computed as where f i is the according vegetation fraction.
Equation (3a) has been derived for plants that cover a significant cross section of the water column that allows for a constant velocity field within the vegetation. This requirement allows for a straightforward separation into one part of the flow parameterized by the vegetation and the logarithmic profile above the

The Dynamic Vegetation Model
Our dynamic vegetation model consists of different vegetation modules divided into colonization, growth, and mortality that parameterize a dynamic, that is, temporally and spatially varying, generic salt marsh species by several literature-based vegetation rules (Figure 3). The interacting processes cause a Note. The parameters t 1 , t 2 , and t 3 are onset of growth season, onset of maximum biomass in summer, and onset of reduced biomass in winter, respectively. a Bouma et al. (2013). b Cooper (1982). c Davy et al. (2001). d Nehring and Hesse (2008). e Poppema et al. (2017). temporal-physical plant variation by growth and temporal-spatial dynamics induced by mortality. The dynamic vegetation model is based on the riparian vegetation model of Van Oorschot et al. (2017), but the intraannual salt marsh growth and treatment of the periodic tides is novel, as previous studies simplified the hydromorphodynamic stresses by the periodic tides Lokhorst et al., 2018).

Colonization, Growth, and Aging
Seedling establishment occurs at the second ecological time step of each ecological year. The spatial distribution of the seedlings is based on the mean of the computed water levels in the HM-model for the precedent ecological time step. A cell was defined as flooded when the water depth within one ecological time step surpassed 0.02 m to account for soil saturation and evapotranspiration (Hughes et al., 2001;Miller & Zedler, 2003). The cells that were flooded and subsequently dried again are determined as locations where seeds establish. This method is based on the assumption that seeds are dispersed by the tide, which is the main reproduction strategy of fast-colonizing pioneer marsh vegetation such as Spartina anglica and Salicornia ssp. (Schwarz et al., 2018). The vegetation is assigned to the cells by the means of an initial fraction of a max- Note. We investigated the ecomorphodynamic model (V) for three bathymetries (2000,2006,2011) and reference scenarios without vegetation (R) for the same bathymetries. The V2000 scenario was also run for the entire period of 2000-2015 (V2000L) to assess the long-term effects of morphological change and vegetation, complemented by a reference model run without vegetation (R2000L). Static vegetation maps (S) were computed to isolate eco-engineering effects.
Model runs for two coarser grids on the bar (G1Bar, G2Bar) and spanning the entire estuary (G2Est) were conducted to investigate upscaling and grid size effects. imum density for pioneer salt marsh that was typically found in literature (Table 1). We assume a generic perennial pioneer species, based on Spartina anglica and Salicornia ssp., whose underground biomass survives in winter and regrows in the subsequent year. Tussock and patch formation smaller than the grid resolution was ignored here. The physical plant properties of the modeled salt marsh vegetation are determined by its growth function which incorporates seasonal plant dynamics (Ibañez et al., 1999). Seasonal dynamics are expressed in the variables that change within one or several years of plant growth, that is, shoot height, root length, and stem diameter (see Table 1). The vegetation settles in the beginning of the growth season (at time t 1 ) and grows linearly with time to a maximum height and stem diameter (at time t 2 ). Between t 2 until the end of the growth season (at time t 3 ) both remain constant. Belowground biomass is calculated as linear growth of the root between t 1 and t 3 until the maximum root length. With the end of the growth season (at time t 3 ) we assume a reduction in stem height to include decay and standing dead biomass of the plant. Root length and diameter remain constant. The patches surviving the first year age into a different life-stage, from seedling to mature vegetation. For mature vegetation a different growth function is defined that allows them to grow larger and defines a larger bulk drag coefficient to account for leaves. The existence of multiple life-stages in one single grid cell is possible, where the maximum sum of all vegetation fractions in one cell equals one. A fraction represents the relative space in a covered grid cell. The generic salt marsh species is defined for all life-stages by the values in Table 1.

Mortality
While the growth rules define the temporally changing physical plant properties in their initial locations (Figure 3  erosion and sedimentation (Figure 3). For the hydrodynamic stressors (i and ii) we use a linear dose-effect relation where the mortality rate increases with increasing pressure: m is the slope of the linear function, x the stressor strength, and b the intercept. Linear mortality was chosen as it represents the plant's resilience to beginning pressures (Holling, 1973). The inundation thresholds are the same for both life-stages, but velocity resilience increases with age (Cao et al., 2018). We chose the 95 th percentile of the maximum flow velocity of each tidal cycle as the representative velocity for plant mortality to avoid unrealistically high velocities at very small water depths due to numerical effects. Thresholds and derivations were taken from literature and estimated based on empirical studies (see Table 2).
In contrast to a gradual mortality rate for hydrodynamic pressures, the morphological pressures of burial and scour are assumed to be fatal to the entire plant fraction if the critical length l crit of the root or the shoot is exceeded by the erosion or sedimentation rate, respectively (equations (7a) and (7b)). The critical root length is defined as 10% of the root, while the critical shoot length is 100% of the shoot. As root length is growing with time, the resilience of the plant against scour increases with age, but resilience against burial changes with the seasonal variations of the shoot length. Mortality is calculated as follows: for which Table 2 summarizes the mortality parameterization.

Scenarios
We performed four steps (Table 3). For all scenarios intertidal area is defined as bed elevations above mean water.
First, we analyzed the performance of the ecomorphodynamic model nearly independently of the quality of long-term morphological modeling. To this end, the ecomorphodynamic model was run for 4 years on three different initial bathymetries derived from data in 2000, 2006, and 2011. These models are further referred to as vegetation scenarios or V-scenarios compared to reference scenarios without vegetation, further referred to as R-scenarios. These specific years are compared to the available ecotope maps in 2004, 2010, and 2015, which span the development from sparse to dense salt marsh on the tidal bar. We kept our water level boundaries the same for all three runs and their reference runs. We carried out a cell to cell comparison between the ecomorphodynamic model results and the ecotope maps to quantify performance of the vegetation predictions using the MATLAB software package (version 2016a). We binarized the vegetation maps and classified correspondence into four categories: correctly predicted salt marsh, predicted but not observed (false positives), not predicted yet observed (false negatives), and correctly predicted absence of salt marsh. These categories were each summed and divided by total intertidal area. Second, we ran the ecomorphodynamic model over the entire 15 year period and compared it to a morphodynamic reference run without vegetation (further referred to as long term or L-scenario). We analyzed vegetation density and distribution along bed elevation and compared the results to predictions by a static model approach (S-scenarios), where vegetation presence is directly predicted from the inundation time thresholds used in the ecomorphodynamic model.
Third, to control for eco-engineering effects, we plot differences maps of the main parameters that control vegetation growth, hydroperiod, 95 th percentiles of flow velocity and bed level, for the two bathymetries of 2000 and 2011. These maps are derived by subtracting the reference scenario (R) from the scenario with  vegetation (V), hence they show vegetation effects. Here, the last tidal cycle is compared, which means that bed level changes are a cumulative result throughout the entire simulation and represent a general trend, whereas the hydrodynamics are the end result and can directly be linked to the present vegetation cover.
Fourth, we applied the ecomorphodynamic model on the entire estuary scale on the coarser grid to gain insights into the large-scale applicability and grid size dependence. We specifically compare vegetation cover and location, the vegetation age distribution and quantify distributions of bed level, and levels at which salt marsh formed as hypsometric curves. On the other hand, the V2000L scenario over 15 years (Figure 4e), shows dense salt marsh cover at the center of the intertidal bar with a much smaller spatial extent than the 4 year simulations and the ecotope maps (Figures 4a-4d). The only difference with the short-duration scenarios is the long-term development of morphology, indicating that the HM-model did not elevate the bar as much as evidenced by the observations. Two main differences between model and data emerge: all model simulations predict salt marsh to settle in the western part of the bar, which is not found in the ecotope maps (Figures 4a-4e and 5). Dense vegetation observed in the southeastern part of the bar is not completely captured by the model.

Vegetation Establishment in Response to Hydromorphological Conditions
To quantify correspondence and mismatch between predicted and observed salt marsh development we imaged the cell to cell comparison as relative binarized maps ( Figure 5). The ecomorphodynamic model predicted around 60% correctly in all three V-scenarios. The total vegetated area in the model was larger than that in the ecotope maps, mainly due to the (false positive) prediction of salt marsh on the west of the bar. However, the vegetation cover was predicted increasingly well with the newer bathymetries in the V2006 and V2011 scenario. Most of the salt marsh in the middle and on the eastern side of the bar was correctly modeled, while the model underestimated vegetation cover at the marsh edge. To understand how the observed vegetation patterns linked to bathymetry of the respective scenario, we look at the relative bed elevation distribution of the V-scenarios compared to the R-scenarios and the corresponding vegetation cover ( Figure 6). The bed elevations of the bar increase over time from the V2000 to the V2011 scenario as shown after normalization to mean water level ( Figure 6). The bed elevation increase causes a shift of the distribution mode from around 0.55 rMW (relative to mean water level) for V2000 to 0.65 rMW for V2006 and 0.7 rMW for V2011, respectively. Moreover, the maximum elevations of the bar are increasing over time while the area with lower elevations reduces around the lower limit of vegetation occurrence. On the other hand, the distribution of bed elevations of the long-term run V2000L of 15 years is rather similar to the V2000 scenario but much different from the bathymetry in the R2011 scenario that it targeted ( Figure 6). This is the likely explanation for the observation that the vegetation did not extend as much in this model run as it did in the shorter duration vegetation scenarios with the higher observed bed levels.
The development of the average vegetation cover along relative bed elevation for the dynamic V-scenarios (Figure 7) shows that the increase over time is mainly due to the increase in higher elevated bar areas. For Figure 11. Grid size dependency of inundation period, peak velocity (95%), and vegetation settlement for three different grid sizes and bathymetry of V2011 (rows). The increase in grid size increases the inundation period, affecting vegetation cover (columns). Velocity is not significantly affected by the change in grid size. We compare the hydrodynamics and vegetation cover for two different flooding-drying threshold (FD threshold) for the vegetation calculations along the different grid sizes, showing that the increased grid size can be compensated by an increased FD threshold (lower two rows). all bathymetries the area occupied by vegetation increased gradually with elevation, of which the maximum cover is found at the highest elevations. In comparison, the S-scenarios predicted vegetation presence from upper and lower thresholds in inundation period. The predictions are in general agreement with the V-scenarios for the range of elevations that became vegetated, but the static scenarios could not predict the trend in vegetation cover with elevation. As the vegetation cover determines the resulting hydraulic resistance we look into their distributions below (Figure 7).
Modeled cover is strongly related to life-stage. The majority of juvenile vegetation is located between 0.45 and 0.75 rMW and mature vegetation between 0.55 and 0.85 rMW (Figure 8a). This distribution is directly linked with the density classes in Figure 4b, showing that the sparse cover is mostly made up of seedlings while the dense cover occurs mainly on high-lying cells with an established mature salt marsh. Roughly two thirds of the surface of colonized relative elevations (between 0.55 and 0.75 rMW) is covered both by sparse and dense and juvenile and adult vegetation. This indicates that cells with mature vegetation are subject to mortality and can partly be recolonized by seedlings. Only at the eastern tip, where cells were only covered by mature plants, we found a 100% coverage and cells could not be recolonized (Figure 8b).
We found a relation between cover density and mortality causes through the ecological/morphological year (Figure 9). For the V2000 scenario, at each grid cell the mortality causes were extracted and added up to calculate the relative contribution to the entire mortality at each time step over one ecological year. If several mortality causes provoked dying of the plant patch, they were classified as combined pressures. The total cover increased at the beginning of the year, which can be attributed to colonization by new seedlings resulting in sparsely covered cells. Over the course of the simulation parts of the newly established vegetation died, due to high inundation rates, while the mature dense vegetation remained constant. The strongest mortality occurred within the first time steps after colonization as the seedlings were located at exposed locations where all pressures were exceedingly high. This process leads to a dynamic equilibrium of colonization and mortality of the vegetation with time. Interestingly, for higher bed elevations (V2006 and V2011) the ratio dense/sparse is increasing, while in V2000 it is around one (indicated by the similar locations of the dashed and dotted line, Figure 9). This shows that the stability of the marsh is very low if the bed elevation is not sufficiently high. In comparison, the static vegetation scenario does not account for mortality; hence the cover is constant with time.

Eco-Engineering Effects of the Dynamic Vegetation
While hydromorphodynamic conditions influence vegetation cover and density, the vegetation in turn influences these conditions. These eco-engineering effects on hydromorphodynamics were calculated as differences with the reference runs (R-scenarios) without vegetation ( Figure 10).
Three contrasting effects emerged. As expected, the flow velocity magnitude on the bar is reduced through vegetation cover. Second, the morphological difference is fairly limited, with the largest effect occurring in the V2000 scenario when the salt marsh had the smallest extent and nearly absent accretion on the marsh. The third and unexpected result is that the inundation duration increases with salt marsh extent, meaning that the presence of vegetation increases the inundation stress for the vegetation. As the bathymetry of the bar increased in elevation and vegetation cover, all effects became more pronounced (Figure 10).
Flow velocity in and around the vegetated area is considerably lower than without vegetation, up to 0.4 m/s. On the other hand, flow velocity is higher in deeper channels developed in the presence of vegetation, especially in V2011 where the cover was largest.
The indirect feedbacks between vegetation and morphology are limited (Figure 10), except at the northern tip of the bar close to the flood channel where the flow is deviated during flood. This causes high dynamics in the sedimentation and erosion pattern. The flow is accelerated along the marsh edge which leads to erosion next to the marsh and immediate sedimentation next to eroded cells. This effect is especially pronounced on the lowest bathymetry (V2000), which is more regularly flooded and a change in hydrodynamics by vegetation can have a potentially larger effect. Bed level change is fairly limited considering the large vegetation cover on the bar.
The larger roughness due to vegetation led to a longer residence time of the water on the bar and thus increasing the hydroperiod. This effect was increasingly pronounced with higher bar elevations due to smaller relative water depth and increased vegetation cover, reducing the drainage after high water. This is the opposite of the idea that the plant is changing its environment to its favor. The increase in inundation was also observed away from the marsh indicating a larger regional effect of the marsh.

Numerical Effects and Test of Upscaling
The grid size resolution is an important aspect for the simulation of vegetation cover. An analysis of two coarser grids (runs G1Bar and G2Bar) shows that the vegetation cover reduces with grid size; hence, there is a relation between numerical effects and upscaling to larger, coarser grid cells. The decrease in vegetation cover is explained by the increase in inundation period for the coarser grid cells (Figure 11). The velocities were however decently distributed over the different grid sizes. To compensate for the increased inundation period we raised the flooding-drying threshold from 2 to 10 cm for the vegetation calculations ( Figure 11). This led to a realistic vegetation cover on the coarsest grid for the bar (Figure 11, coarse, FD-threshold:0.1 m).
To test the possibility for upscaling the vegetation model to a coarse grid, we applied a vegetative flooding-drying threshold of 10 cm (G2Est scenario). The results in comparison with the ecotope map are adequate ( Figure 12). The salt marshes along the estuary margins were well reproduced, as were the large bars "Hooge Platen" and "Walsoorden." The salt marsh of "The drowned land of Saeftinghe" was underrepresented with a lower cover than in the ecotope maps, while the ecomorphodynamic model falsely predicted plant cover on some smaller bars. The areas with false cover predictions are classified as low-dynamic high and supralittoral or muddy littoral in the original ecotope map.

Discussion
The vegetation cover acquired with our novel ecomorphodynamic model compares well with the ecotope maps both for salt marsh location and density. Literature-based values of vegetation properties are sufficient to predict the locations and cover of the generic salt marsh with an accuracy of 60% after a simulation period of four morphological years. Hence, a calibration of the vegetation parameters in the dynamic vegetation model for different bathymetries is not necessary, which makes the model unique and opens up possibilities to evolve a generic code for a variety of systems in the future.
We found salt marsh density and age (expressed as life-stages) to positively correlate (Figure 8), which agrees with previous observations in salt marsh zonation (Adam, 1993). In contrast with literature, our model results show an increase of the hydroperiod next to the established salt marsh, suggesting an unexpected negative eco-engineering effect ( Figure 10). More specifically through the salt marsh-induced increase in hydroperiod, salt marsh expansion becomes limited adding a new dimension to scale-dependent vegetation-landform interactions. Possibly, this opens up niches for other species not yet included. Finally, we show that the grid size affects the calculation of the drying and flooding hydrodynamics and thus the prediction of vegetation patterns.

Density Gradient-Plant Age and Fast-Slow Colonizers
Our model results show a vegetation density gradient along the bed elevation (Figure 7). This is a typical pattern observed in pioneer marshes, where the high marsh is characterized by dense and mature vegetation while limited seedling survival and lateral expansion cause sparse cover at low elevations (Gray & Bunce, 1972). The patterns in our model emerged from the dynamic, literature-based rules that are prescribed for the two life-stages, where seedlings are smaller and more susceptible to velocities. The young plants need a disturbance-free period to grow sufficient roots that can withstand the flow (Balke et al., 2012;Wiehe, 1935). As soon as sufficient time has passed, the vegetation grows and develops a strong root system that prevents erosion, which is parameterized by the mature life-stage (Cao et al., 2018). Our model results indicate that the main eco-engineering effect for plant survival is that the increase in plant height and stem width reduces the flow strength, and at the same time increases plant resilience. This reduced flow also causes low erosion rates on the bar, which renders possible stabilizing effects of plant roots unimportant ( Figure 10). Consequently, the altered environmental conditions by these sparse seedlings facilitate the survival of the next generation of seedlings.
At the same time, since the ecomorphodynamic model allows different life-stages within the same grid cell it also incorporated shielding of younger plants by the older and larger plants . On the other hand, the life-stage-specific stress tolerance of the seedlings also causes increased mortality at these lower elevations. which results in the simulated density gradient as observed in field studies (Christiansen et al., 2000;Van der Wal et al., 2008).This effect can also be observed at small-scale salt marsh-mudflat features, such as unvegetated tidal channels caused by flow concentration between vegetation patches Van Wesenbeeck et al., 2008). The degree of flow concentration and channel emergence is a function of the flow field, vegetation density, and water depth. Our model reproduces flow concentration through the interplay between eco-engineering effects and life-stage-dependent growth and mortality ( Figure 10).
Salt marshes are characterized by various primary colonizers (Schwarz et al., 2018), which we only indirectly modeled. First, species establish that have low densities and short life spans, which reduce the flow and enhance bed accretion. These fast colonizers facilitate establishment of successional species by their eco-engineering activity. Slow colonizers usually have higher densities and are perennial, meaning that they start regrowing in the subsequent year. Throughout their life span both strategists can cooccur at similar elevations on the intertidal flat (Suchrow & Jensen, 2010). We indirectly modeled these traits with our generic salt marsh species as a mix of fast and slow colonizers through specification of the life-stages. In doing so, we first allow an establishment of the rather small fast colonizers on the bare sandflat with sparse densities. In locations with sufficient shielding through elevation or eco-engineering the plants survive their first year and become more resilient plants that can grow higher and in larger densities. At the same time, they become more resistant to velocity stress, which is a main trait of slow-colonizing plants. This means that the modeled generic species represents the properties of different pioneer species as Salicornia ssp. and Spartina anglica through the life-stages.

Feedback Between Eco-Engineering Effects and Salt Marsh Vegetation
Hydroperiod is the main control of vegetation growth and sedimentation, which depends equally on the water levels and the bed elevations of the system (Adam, 1993;Allen, 2000;Kirwan et al., 2010;Mudd et al., 2004;Reed, 1990). However, water levels and bed elevations are co-dependent, vegetation reduces flow velocity in dense marshes, which influences hydroperiod and leads to higher sedimentation rates at the marsh edge rather than in the inner marsh (Townend et al., 2011). Also, Neumeier and Ciavola (2004) observed that canopies act more as erosion protection than enhancing accretion under normal conditions. This agrees with our results, where bed level changes mainly occur adjacent to the marsh, while the short hydroperiod and the low sediment mobility of the single grain size within the marsh caused insignificant sediment transport.
In contrast with previous studies, a new negative eco-engineering effect emerged. While locally the salt marsh reduced the velocity magnitudes, making conditions more favorable for further settling, the inundation period was increased by the presence of extensive vegetation cover, which acts as a constraint on its spatial extent ( Figure 10). Especially at the lower elevations away from the dense marsh edge, this effect causes mortality of young vegetation. Consequently, this new eco-engineering effect directly limits the marsh's growth and hampers salt marsh expansion in sandy dynamic systems on marsh scale. It was already

Journal of Geophysical Research: Earth Surface
10.1029/2019JF005092 stated by Nepf (1999) that effects on flow increase with distance from the marsh, who were further supported by the theory of large-and small-scale-dependent feedbacks between vegetation and mudflat (e.g., Schwarz et al., 2015Schwarz et al., , 2018Van Wesenbeeck et al., 2008). Similarly, we found eco-engineering effects on flow rates at the marsh scale. Interestingly, even though the cover increased gradually with elevation in the ecomorphodynamic model, the positive eco-engineering effect on the velocity magnitudes is relatively strong (Figure 10 : 2004, 2010, 2015). Contrarily, the "marsh-scale negative eco-engineering effect" only seems to become visible after a certain salt marsh size developed (Figure 10: 2010, 2015).
To quantify the importance of the eco-engineering effects on vegetation patterns we compared the results of our ecomorphodynamic model to a static model for the lowest bed elevations in 2000 and highest bathymetry of 2011 ( Figure 13). The static model prescribes vegetation growth based on thresholds for inundation period. We determined the mean inundation period over the 4 years of a run without vegetation cover and prescribed vegetation based on the benchmarks found in literature. We see that for the same threshold of 0.3 inundation period the static vegetation cover is larger for both S-scenarios than the corresponding ecomorphodynamic model results, which can be directly attributed to the marsh-scale negative eco-engineering effect of the inundation period. To quantify their difference we add a lower threshold of 0.2 to the static maps which produces a similar pattern as the ecomorphodynamic model with 0.3. The marsh-scale negative eco-engineering effect already shown in Figure 10 led to higher mortalities that reduced the vegetation cover in the ecomorphodynamic model compared with the static model, limiting the predictive capabilities of static model approaches. At the same time, positive engineering effect for velocities were found where dense vegetation can expand further on the bar while it is shielded from present plants (corresponding with Figure 10). This supports our theory that velocity reduction strongly facilitates the survival of seedlings and allows for marsh expansion, even without large sediment supply.
The results stress the importance of a dynamic vegetation model, since vegetation alters the pressures causing its survival. While some factors are improved by the plants, the most important factor, inundation period, is not positively eco-engineered in the absence of fine sediment trapping and significant organic material production. As a result, vegetation predictions require consideration of dynamic feedbacks between newly established vegetation and environmental factors, including a detailed representation of their physical appearance that controls the magnitude of their eco-engineering effect. Static predictions do not capture the variations in eco-engineering effects. While being unsuitable for predictions of new salt marshes, static predictions can give insights into the impact of vegetation on morphology in hindcasting scenarios. This finding is especially important for modeling studies predicting future colonization of salt marsh, quantification of coastal protection, and management strategies on the system scale and along the fluvial-tidal transition.

Effects of Limitations in the Morphological Model
A comparison with a continuous run simulating morphologic development over 15 years showed a bigger disparity between the model and the ecotope maps ( Figure 4). We believe that this difference is not caused by shortcomings in the dynamic vegetation model, but is related to lower bed accretion rates than observed in reality, possibly caused by simplifications of the HM-model, that is, the choice of the sediment fraction or disregarding dredging activities. Even though the Western Scheldt is predominantly sandy, natural variations in sediment fractions as well as redistribution of sediments by dredging and dumping possibly lead to larger sedimentation rates than those predicted by the model. The HM-model includes one sand fraction of 200 μm that is representative for the main channels but neglects finer fraction typically observed in more sheltered areas (Braat et al., 2017;Gray & Bunce, 1972), which potentially are important for their vertical buildup. Vegetation uses organic accumulation and fine sediment trapping as two mechanisms to increase bed level elevation (Braat et al., 2017;D'Alpaos, 2011;Le Hir et al., 2007;Turner et al., 2002). Due to the choice of the sediment and the vegetation growth strategy, both processes are not yet accounted for in the ecomorphodynamic model. This can explain the low bed level differences between the R-and V-scenarios ( Figure 7) and the reduced bed accretion rates compared with reality.
The other factor controlling marsh expansion and bed accretion is sediment supply (Friedrichs & Perry, 2001;Mariotti & Fagherazzi, 2010). After a pilot dumping of 500,000 m 3 in 2004, dumping adjacent to the bar has been carried out extensively from 2010 onward with yearly between 500,000 m 3 to 1,300,000 m 3 of dredged material (Leys et al., 2006;Plancke et al., 2017). The artificial sediment supply by the redistribution of the dumped sediment can have had a significant part in the observed increase in bed elevation between At the same time, vegetation covers a larger area on the bar compared to the ecotope maps. This could be linked to the roughness distribution on the bar as the HM-model only considers sand. Mud that can be found on the bar was not included in the calibration due to large computational expenses and prior calibration efforts, but would result in lower roughness values in reality. Sensitivity runs with a lower bottom roughness in both domains resulted in higher velocities on the western side of the bar, which removed newly established vegetation. Consequently, roughness effects induced by sediment type are a possible contributor to reduce vegetation abundance. At the same time, at the northeastern bar margins vegetation cover is lower in the model ( Figure 5). The model is not able to predict the bed forms that are characterizing that area, causing higher flow and mortality in the highly dynamic areas near the flood channel.
Sediment type also controls vegetation density. Van Hulzen et al. (2007) showed that sediment accumulation and vegetation density were enhanced on mudflats compared to sandflats. Parameterization of sediment dependence of vegetation in the dynamic vegetation model can give new insights into their establishment mechanism and offer new possibilities for coastal sediment management. For example, modeling with sand and mud fractions in combination with the dynamic vegetation model could elucidate under what conditions mud flats appear that facilitate plant settling, or plants settle that cause mud to accrete.
The entire estuary model with vegetation largely reproduced the cover in the ecotope maps ( Figure 12). However, the model predicts a smaller vegetation cover at the large Drowned land of Saeftinghe marsh than in the ecotope maps, which is due to a combination of two factors. First, in the original flow model the Saeftinghe area was calibrated with a high constant hydraulic roughness due to the pronounced vegetation cover (Vroom et al., 2015). This led to high hydroperiods as the water was unable to flow off the marsh after flooding, causing high and unrealistic salt marsh mortality in our ecomorphodynamic model. Second, the natural bed elevations of the Saeftinghe marsh are high compared to the bars. This is due to the high age of the marsh which allowed vegetation-induced sedimentation to such high bed elevations that the marsh cannot be flooded every tidal cycle (Stark et al., 2017).
At the same time, some of the bars were predicted to have salt marsh establishment while the observations show bare soil. This can be linked to growth-hampering processes lacking in the ecomorphodynamic model, such as shipping waves, strong secondary currents, and highly three-dimensional flows in the sharpest bend in the middle of the model grid. Interestingly, most of the falsely predicted marsh is linked to the low energetic classes or supralittoral ( Figure 12). These areas have a high bed elevation and are already low in dynamics, hence, potentially suitable for vegetation establishment.
Stochastic (e.g., Schwarz et al., 2015Schwarz et al., , 2018Temmerman et al., 2007) or deterministic model approaches (e.g., D'Alpaos et al., 2012;Kirwan & Murray, 2007) both allow for a good representation of natural salt marsh patterns. Both approaches are able to simulate salt marsh development: while the former is often used to describe emergent properties such as tidal channel formation, the advantage of the latter is to be able to compare governing parameters through deterministic scenarios. Patch formation of salt marsh vegetation is not represented in the current model as we focus on large-scale morphology and salt marsh growth requiring larger grid sizes. In spite of the model neglecting small-scale vegetation and morphological features, we show that realistic vegetation patterns emerge that give new insights into the interactive effects between salt marsh growth and morphology.
Further investigations are needed to link multiple species effects on sediment. We suspect that sediment supply is the limiting factor for bed level accretion in the model, which creates the need of adding fine sediment supply and dumping to the model for future model runs. The effects of other ecological features of vegetation need to be tested to understand the parameters needed to accurately model vegetation effects in fluvial and coastal environments and enhance morphological modeling.

Conclusion
Without calibration, our ecomorphodynamic model correctly predicted the decadal salt marsh settling and expansion trend on the bar of Walsoorden based on vegetation parameters reported in literature.
The ecomorphodynamic model produced diverse vegetation patterns through density gradients and eco-engineering effects. Density of the marsh varies spatiotemporally due to the presence of multiple life-stages that have different physical plant properties and resilience. These life-stages also allow distinctions between fast and slow colonizers, which opens up possibilities to test a variety of ecological concepts and ecogeomorphological interactions.
While hydromorphic conditions determined the vegetation settling and initial mortality patterns, two eco-engineering effects of vegetation emerged.
(1) The reduced flow velocities due to hydraulic resistance allow for an expansion of the salt marsh area.
(2) Extensive vegetation cover increased the inundation period causing higher mortalities, which resulted in a negative eco-engineering effect. As a result, in sand-dominated systems salt marshes emerge from an equilibrium between positive and negative eco-engineering effects. A static vegetation presence model based on hydrodynamic thresholds mispredicted the observed vegetation patterns because it lacked the feedbacks between vegetation and hydromorphodynamics. Consequently, spatial pattern and stability of the marsh depend on a combination of shielding by vegetation as well as plant resilience, which depends on the age of the marsh.
When predicting vegetation cover along the entire Western Scheldt estuary, grid size effects need to be accounted for to accurately represent the hydromorphodynamics. The underrepresentation of the drying and flooding process on the coarser grids can be compensated by a different water depth threshold for inundation.
Our results support the hypothesis that salt marshes in estuaries in fact control not only local processes but also the marsh scale. This suggests that vegetation affects the large-scale pattern of channels, bars, and overall planform of estuaries.