Biogeomorphic modeling to assess the resilience of tidal-marsh restoration to sea level rise and sediment supply

Abstract. There is an increasing demand for the creation and restoration of tidal marshes around the world, as they provide highly valued ecosystem services. Yet restored tidal marshes are strongly vulnerable to factors such as sea level rise and declining sediment supply. How fast the restored ecosystem
develops, how resilient it is to sea level rise, and how this can be steered by restoration design are key questions that are typically challenging to assess due to the complex biogeomorphic feedback processes involved. In this paper, we apply a biogeomorphic model to a specific tidal-marsh restoration project planned by dike breaching. Our modeling approach integrates tidal hydrodynamics, sediment transport, and vegetation dynamics, accounting for relevant fine-scale flow–vegetation interactions (less than 1 m2) and their impact on vegetation and landform development at the landscape scale (several km2) and in the long term (several decades). Our model performance is positively evaluated against observations of vegetation and geomorphic development in adjacent tidal marshes. Model scenarios demonstrate that the restored tidal marsh can keep pace with realistic rates of sea level rise and that its resilience is more sensitive to the availability of suspended sediments than to the rate of sea level rise. We further demonstrate that restoration design options can steer marsh resilience, as they affect the rates and spatial patterns of biogeomorphic development. By varying the width of two dike breaches, which serve as tidal inlets to the restored marsh, we show that a larger difference in the width of the two inlets leads to higher biogeomorphic diversity in restored habitats. This study showcases that biogeomorphic modeling can support management choices in restoration design to optimize tidal-marsh development towards sustainable restoration goals.


: Vegetation dynamics model scenarios (i.e., reference vegetation dynamics and two variants, respectively without vegetation and with instantaneous colonization - Table S1). Evolution of the mean platform elevation with respect to the mean high-water level (MHWL) (a-b) and development of the vegetation cover (c-d) in the Northern (a, c) and Southern basins (b, d).  Table S1). Bed elevation 50 years after de-embankment.
The dashed lines delineate the old marsh, the Northern basin, and the Southern basin. The ellipses emphasize a preexcavated channel that has disappeared (a-b) or survived (c), depending on the vegetation dynamics. All figures are rotated by 43° clockwise, as compared to Fig. 2c. Figure S3: Vegetation input parameter model scenarios (i.e., reference vegetation dynamics and four variants, respectively with low and high establishment probability (a, c), and with low and high lateral expansion rate (b, d) - Table S1). Evolution of the mean platform elevation with respect to the mean high-water level (MHWL) (a-b) and development of the vegetation cover (c-d).     Figure S1 to Figure  Instantaneous colonization Table S5.

S1 Biogeomorphic model
We have developed the biogeomorphic modeling framework Demeter to simulate explicitly the feedbacks between hydrodynamics, morphodynamics (Sect. S1.1) and vegetation dynamics (Sect. S1.2). This is a multiscale approach, in which the vegetation dynamics is computed at much finer resolution than the hydro-morphodynamics (Fig. 1), requiring the development of specific multiscale coupling techniques to preserve subgrid-scale heterogeneity while information is exchanged between the hydro-morphodynamic and vegetation modules (Sect. S1.3 and S1.4). The specific setup for our study site is detailed in Sect. S1.5.
Telemac-2d solves the depth-averaged shallow water equations in a two-dimensional horizontal framework (Hervouet, 2007) to simulate fluctuations of the water depth ℎ and the depth-averaged flow velocity : ℎ + ⋅ (ℎ ) = 0 (S1) where is the time, is the spatial differential operator, is the gravitational acceleration, is the water surface elevation above the reference level (NAP), is the diffusion coefficient, is the bed shear stress, is the vegetation resistance force per unit horizontal area, and is the water density. The bed shear stress is computed with the Manning formula: where the Manning coefficient is empirically derived and depends mainly on bed roughness.
The vegetation resistance force is modeled as the drag force on a random or staggered array of rigid cylinders with uniform properties (Baptist et al., 2007) and depends on the spatial distribution of vegetation provided by the cellular automaton (Sect. S1.4).
Sisyphe solves the depth-averaged advection-diffusion equation to simulate fluctuations of the depth-averaged suspended sediment concentration : where and are the rates of sediment erosion and deposition, respectively. The rate of sediment erosion is computed using the equation of Partheniades (1965): where is the Partheniades constant and is the critical bed shear stress for sediment erosion. The rate of sediment deposition is computed using the equation of Einstein and Krone (1962): where is the sediment settling velocity. The existence of a threshold shear stress below which sediments remain in suspension is debated in the literature. Here we follow one of the well-established arguments that such threshold does not exist, and that it rather represents a threshold for erosion of freshly deposited sediments (Winterwerp, 2007). This approach agrees with field observations in the Chesapeake Bay (Sanford and Halka, 1993) and is often adopted in recent biogeomorphic models (e.g., Adams et al., 2016;Bryan et al., 2017;Mariotti, 2018;Zhang et al., 2019;Brückner et al., 2020).
The evolution of the bed is computed as follows: where is the bed surface elevation above the reference level (NAP), is the morphological acceleration factor (Sect. 2.1) and is the sediment dry bulk density. The bed is composed of two layers: the fresh layer at the surface and the compacted layer underneath. Their evolution obeys the following rules: (i) each layer is characterized by different values of and , (ii) erosion of the compacted layer only occurs where and when the fresh layer is locally empty, (iii) deposition only occurs on the fresh layer, and (iv) there is no sediment flux between the two layers.

S1.2 Cellular automaton (vegetation dynamics)
As vegetation module, we use the cellular automaton implemented in Demeter. A cellular automaton consists of a regular grid of cells, each one with a finite number of states (here, either bare or one of the considered vegetation species). Cells can change their state in discrete time steps, depending on their neighborhood state and a set of simple stochastic transition rules (Balzter et al., 1998).
Our cellular automaton is implemented as a hierarchical model, where higher-rank species are stronger competitors able to outcompete lower-rank species. In our model, higher-rank species can displace lower-rank species, but not the other way around. Lower-rank species can only colonize after higher-rank species have died off. On the long term, high-rank species will therefore always outcompete lower-rank species in their own niche.

S1.2.1 Establishment
Establishment is the transition from bare state 0 to any vegetated state . The probability of establishment for species is evaluated as: where is the background probability of establishment for species , and are stress functions of the environmental variables (Sect. S1.2.5).

S1.2.2 Succession
Succession is the transition from any vegetated state to another vegetated state > (e.g., from pioneer to climax vegetation). The probability of succession , from species to is evaluated as: where , is the background probability of succession from species to .

S1.2.3 Stress-related die-off
Stress-related die-off (or simply die-off) is the transition from any vegetated state to bare state 0 due to environmental stress. The probability of die-off for species is evaluated as follows:

S1.2.4 Annual die-off
Annual die-off is the transition from any vegetated state to bare state 0 due to the natural cycle of annual species. The probability of annual die-off for species is evaluated as follows: where is the background probability of annual die-off for species .

S1.2.5 Stress functions
Stress functions (Sect. S1.2.1 to S1.2.3) can be of two shapes. When vegetation is only affected at high (resp. low) values of an environmental stressor, and not below (resp. above) a certain threshold, we use the Hill function, which varies from 0 to 1 following: where is the environmental variable, is the threshold around which the transition from 0 to 1 occurs, and is a parameter that controls the shape of the function. The function decreases from 1 to 0 if < 0 and increases from 0 to 1 if > 0. The transition from 0 to 1 becomes steeper for increasing | |.
When the range of optimal conditions is confined between a low and a high threshold value, we use the Brière function: where 0 and 1 are the low and high thresholds, respectively, and is a coefficient used to rescale the function, so that its maximum value is 1: The different environmental variables used for the stress functions are the hydroperiod, the bed elevation gain and loss, and the binned shear stress (Sect. S1.3.1).

S1.2.6 Lateral expansion
Lateral expansion is the transition from any state (bare or vegetated) to any vegetated state > resulting from the presence of at least one neighboring cell of state . The recruitment process is here quite different than for the other processes. It is defined by the mean expansion rate , which determines the number of iterations of the cellular automaton. For each iteration, the probability of recruitment by lateral expansion is where Δ is the grid resolution of the cellular automaton. With this stochastic approach, even though the mean expansion rate is constant, the actual expansion rate varies in space and time. The number of iterations is determined so that > + 2 2 (S17) where the maximum expansion rate and the variance of the expansion rate 2 are calculated as follows: As each species can have a different mean expansion rate, and hence a different number of iterations, we use the highest number of iterations among all species.

S1.2.7 Computational sequence
The different transition rules of the cellular automaton are scheduled as follows: 1. Annual die-off is applied for each annual species in one single iteration.
2. Establishment, succession, and lateral expansion are applied for all species in an iterative process. The number of iterations is determined based on the mean expansion rates (Sect. S1.2.6). For each iteration, the probabilities of establishment, succession and lateral expansion are rescaled as follows: where is the number of neighboring cells vegetated with the same species at the previous iteration. We use a factor ¼ in Eq. S22, so that the rescaling factor 4 is 1 on average.
3. Stress-related die-off is then applied in one single iteration. S1.3 Coupling Telemac to cellular automaton S1.

Environmental variables
The hydroperiod is the percentage of time during which a Telemac grid node is flooded (i.e., the water depth higher than 0.1 m) between two cellular automaton calls. It varies between 0 (never flooded) and 1 (always flooded).
The bed elevation change Δ is the difference between the final and initial bed elevations between two cellular automaton calls. The bed elevation gain Δ + and the bed elevation loss Δ − are calculated as: The binned shear stress is calculated by classifying flow directions into 8 directional bins (45° each) occurring between two cellular automaton calls. The relative binned time , the binned shear stress ̅ , and the binned water depth ℎ ̅ are respectively the percentage of time, the mean bed shear stress, and the mean water depth when the flow is oriented in the th bin. As bed shear stress and flow directions are especially relevant above certain thresholds of the water depth and the bed shear stress, these binned variables only account for situations when the water depth is higher than 0.1 m and the bed shear stress is higher than 0.1 N m -2 .
The mean water depth ℎ ̅ between two cellular automaton calls is calculated for situations when the water depth is higher than 0.1 m.

S1.3.2 Spatial refinement
We use a linear interpolation to spatially refine the hydroperiod, and the bed elevation gain and loss from the Telemac grid to the cellular automaton grid.
We use the concepts of Voronoi neighborhood to spatially refine the relative binned time and the binned water depth. Each cellular automaton grid cell is associated with its closest Telemac grid node. The Voronoi neighborhood of a Telemac grid node is the ensemble of all associated cellular automaton grid cells. Here, the relative binned time and the binned water depth of a Telemac grid node are passed to all cellular automaton grid cells of its Voronoi neighborhood.
For the binned shear stress, we use a convolution method that allows to account for interactions between flow and subgrid-scale vegetation patterns (Gourgue et al., 2021).
Practically, we first calculate the binned velocity ̅ on the Telemac grid as follows: Then, we use a convolution method (Gourgue et al., 2021) to spatially refine the mean binned velocity from the Telemac grid to the cellular automaton grid. Finally, we calculate the binned shear stress on the Telemac grid as follows:

S1.3.3 Stress function of the binned shear stress
A stress function of the binned shear stress (typically using the Hill function) requires a specific treatment to combine all its components. It is calculated as follows: where is the transmittance coefficient (Sec. S1.4.2), and , , , and are respectively the bulk drag coefficient (Baptist et al., 2007), the vegetation cover (Sec. S1.4.2), the stem density, the stem diameter and the stem height of species .

S1.4.2 Spatial coarsening
The vegetation cover of the species is the percentage of cellular automaton cells of state within the Voronoi neighborhood of a Telemac grid node (Sec. S1.3.2). It varies between 0 (not covered by species ) and 1 (fully covered by species ). The sum of all vegetation covers also varies between 0 (bare) and 1 (fully covered by vegetation).
The transmittance coefficient accounts for the spatial heterogeneity of the vegetation distribution at the subgrid scale (i.e., within a Voronoi neighborhood). In general, hydrodynamic models assume a uniform spatial distribution at the subgrid scale (here, = 1), which leads to considerable overestimation of the flow resistance if the vegetation presents clustered patterns at the subgrid scale (Gourgue el al, 2019). The method to compute the transmittance coefficient builds on the similarity between the Chézy formula in fluid dynamics and Ohm's law in electricity. Taking the analogy further, we recalculate the coarse-scale hydraulic roughness just as the total resistance of an electronic circuit that combines resistors (equivalent to cellular automaton cells in our analogy) connected in series based on previous studies in the same restored tidal marsh area (Maximova et al., 2014;Zhou et al., 2016), the Scheldt Estuary (van Leussen, 1999;Van de Broek et al., 2018) and other intertidal environments (D'Alpaos et al., 2021). They are summarized in Table S3. The suspended sediment concentration at the open boundary and the rate of sea level rise vary according to model scenarios (Table 1).

S1.5.2 Vegetation module
The study site is in the oligohaline zone (0.5 -5 PSU) where Aster tripolium is often observed as the pioneer species, and Scirpus maritimus and Phragmites australis in the marsh interior (Van Braeckel et al., 2008). Their expected encroachment in our study site is further supported by the results of transplantation experiments carried out in nearby tidal marshes.
Aster tripolium is an annual species, which can be found as lower pioneer in calm areas and along creek edges. It colonizes the tidal flats and creek levees every year from seeds, as randomly scattered high density clusters on tidal flats. Although it is regarded as an annual species, part of the established plants can survive and develop for another year. Scirpus maritimus is the dominant perennial species from the low pioneer zone into the middle marsh zone. It is even the only species present in the pioneer zone in several tidal marshes close to the study site. The main mode of colonization on bare tidal flats is via lateral spread of rhizomes (Silinski et al., 2016). Phragmites australis is the dominant species in the high marsh zone. It can form large stands from the high pioneer zone up to the supratidal zone, but it is mostly found above Scirpus maritimus in the middle and high marsh zone. Most seedling establishment occurs within already established vegetation, but very rarely on bare tidal flats, except for the highest areas. Once established, it can often outcompete Scirpus maritimus and colonize vegetated areas by lateral expansion via rhizomes, resulting in clearly visible circular patches within Scirpus maritimus marshes.
The initial vegetation distribution is based on aerial pictures before de-embankment.
Marshes that will be excavated and farmland are considered as unvegetated.
Parameterization of the different stress functions (Sec. S1.2.5) is based on field and flume experiments, remote sensing, literature data and model calibration (Table S4-Table S5).

S2 Sediment accretion on vegetated platforms
Based on digital elevation maps derived from historical topographic surveys in the adjacent marshes of the Drowned Land of Saeftinghe (Fig. 2c) between 1931and 1963(Wang and Temmerman, 2013, we have developed an empirical relationship between mean elevation change on vegetated platforms and mean high-water depth (Vandenbruwaene et al., 2014).
Here, we develop a similar relationship based on model results in the restored tidal marsh, using the same variables over the same time interval (i.e., between years 18 and 50 after deembankment), and we compare it with the empirical relationship derived from observations.
The digital elevation maps derived from historical topographic surveys have a resolution of 20 m. To focus on vegetated platforms and avoid the influence of tidal channels, we only consider vegetated areas that are at least 200 m from tidal channels in the digital maps (Vandenbruwaene et al., 2014). Similarly, as our model results have a resolution of 5 m, we only consider areas that are at least 50 m from tidal channels in the model results.
The Drowned of Saeftinghe is located downstream of the study site, where the sediment input from the Scheldt Estuary is substantially lower. Historical measurements in the period 2001-2012 reveal that the tide-averaged SSC in the estuary is 42 mg l -1 close to the Drowned of Saeftinghe and 63 mg l -1 close the study site (Vandenbruwaene et al., 2014). To account for this 1.5 ratio in sediment input between model and observations, we multiply the observed mean elevation change by 1.5 to obtain the data presented in Fig. S1.

S3 Pioneer vegetation development
We compare our model results with observed rate of spatial expansion of the vegetation cover in the adjacent restored marshes of Paardenschor (Fig. 2c), from the onset of vegetation in 2007 until 2017. We use a series of Google Earth images, and we apply the method of Richardson et al. (2009)

to classify vegetation pixels. Part of the vegetation colonization in
Paardenschor starts from the dikes. Such phenomenon is expected to be of a much lesser influence in our study site. Hedwige-Prosper Polder is about 30 times larger than Paardenschor, hence the average distance to dikes will be much higher. In our analysis, we therefore remove the vegetation development occurring from the dikes.

S4 Channel network characteristics
We compare various geometric properties of the simulated tidal channels with observations in the adjacent marshes of the Drowned Land of Saeftinghe ( Fig. 2c - Vandenbruwaene et al., 2013Vandenbruwaene et al., , 2015. To that end, we have developed a quasi-automatic methodology to extract tidal channel networks and related characteristics from model results. We first identify grid nodes within channels by applying a multi-window median neighborhood analysis (Liu et al., 2015) on the simulated topography, and we compute the unchanneled flow length as the shortest distance to a channel grid node (Tucker et al., 2001). We then retrieve channel edges as multiple polygons by applying the Python function tricontour from the visualization library Matplotlib (Hunter, 2007) on the channel grid nodes. We finally extract the channel network skeleton, defined as the channel centerlines (Fagherazzi et al., 1999), by generating the raw Voronoi diagram of the channel edge polygons (with the Python library Centerline) and applying straightforward threshold rules to simplify it.
We use a virtual topography method to determine the watershed areas along the network skeleton (Vandenbruwaene et al., 2013(Vandenbruwaene et al., , 2015. In terrestrial river networks, watershed areas are exclusively delineated by topographic gradients. For tidal channel networks, however, topographic gradients are small and water flow is mainly determined by water surface gradients (Rinaldo et al., 1999). Alternatively, algorithms designed for terrestrial river networks (here the Python library pysheds) can be applied on a virtual topography built as the sum of the shortest distance to the network skeleton and the distance to the mouth along the network skeleton. For every point along the network skeleton, we can then compute the watershed area and the upstream mainstream length, defined as the longest upstream channel within the corresponding watershed.
Cross-sectional dimensions of tidal channels are traditionally related to the spring tidal prism (D'Alpaos et al., 2010). For tidal marsh channels, however, overmarsh tides that overtop the intertidal platform are more relevant (Vandenbruwaene et al., 2013(Vandenbruwaene et al., , 2015 because maximum channel flow velocities typically occur when the surrounding platform is flooded and drained (Bayliss-Smith et al., 1979;Pethick, 1980;French and Stoddart, 1992).
Here we use the mean overmarsh tidal prism, defined as the mean tidal prism from all overmarsh tides. For every point along the network skeleton, we compute the mean platform elevation of the corresponding watershed. The mean overmarsh tidal prism is then simply the product between the watershed area and the mean overmarsh high-water depth, obtained from all simulated high tides higher than the mean platform elevation.
We generate channel cross-sections along the network skeleton by balancing two constraints: cross-sections must be as perpendicular as possible to the network skeleton and consecutive cross-sections must not intersect each other. Where both constraints can be met, we then compute the channel depth as the difference between the mean channel edge elevation and the lowest cross-section elevation, the channel width as the distance between channel edges, and the cross-section area as the integral of the difference between the mean channel edge elevation and the cross-section elevation.