Intertwined Eco‐Morphodynamic Evolution of Salt Marshes and Emerging Tidal Channel Networks

The formation and development of tidal channels and salt marshes are controlled by complex interactions between hydrodynamics, sediment transport, and vegetation dynamics. Tidal channels affect and, at the same time, are affected by the growth of salt marshes fringing them. The coupled evolution of these morphological units, mediated by vegetation growth, is thus a key ingredient for simulating the behavior of tidal environments. Considering these two factors, we developed a mathematical model to investigate the eco‐morphodynamic evolution of intertidal areas fringing a main channel and of the tidal creeks cutting through them. Model results indicate that vegetation promotes the development of channel networks, leading to more complex channel structures and higher drainage efficiency. Vegetation encroachment influences sediment deposition patterns by trapping sediment in the seaward and middle intertidal areas, while reducing the amount of sediment delivered to landward areas. In the presence of sea level rise, this deficit of sediment enhances the landward‐decreasing trend of the intertidal platform and leads to more isolated vegetation patches. Overall, sea level rise restricts the extension of salt marshes and consequently reduces the effect of vegetation on channel network form and function.

flat when marsh accretion rate is lower than the rate of RSLR, and the marsh platform progressively drowns Kirwan et al., 2010;Marani et al., 2007). In general, the window of opportunity (Balke et al., 2011Hu et al., 2015) for vegetation growth is reduced as the MSL increases and vegetation progressively disappears (Kirwan & Megonigal, 2013;Morris et al., 2002). Salt marshes can also undergo a transition into tidal flats through lateral retreat of marsh edges. These edges are characterized by scarps between the vegetated marsh surface and the bare tidal flat bottom facing the marsh. Wind waves impinge against salt marsh edges, leading to lateral marsh retreat, and erode sediments from the tidal flats, promoting their deepening (e.g., Leonardi & Fagherazzi, 2014;Leonardi et al., 2016;Mariotti & Fagherazzi, 2013;Marani et al., 2011). Deeper and wider tidal flats favor the development of stronger wind waves (Mariotti & Fagherazzi, 2013;Tommasini et al., 2019) and, consequently, enhance wave-induced erosion of the salt marsh border (Bendoni et al., 2016).
The evolution of tidal flats and salt marshes is strongly related to tidal channel dynamics. Both tidal channel and intertidal platform evolution are governed by the interactions between hydrodynamics, sediment transport, morphological changes, and biological effects (Alizad, Hagen, Morris, Bacopoulos, et al., 2016;D'Alpaos et al., 2007;De Swart & Zimmerman, 2009;Gong et al., 2018;Li et al., 2019;Temmerman et al., 2007;Vandenbruwaene et al., 2011;Zhao et al., 2019). The mutual feedbacks (either positive or negative), which characterize these interactions, control the so-called Morphodynamic Loop (Coco et al., 2013;Cowell & Thom, 1994;Wright & Thom, 1977). Given proper initial and boundary conditions, the overall evolution of tidal systems is controlled by the interplay between the various morphological units (tidal channels, tidal flats, and salt marshes). Tidal channels act as preferential drainage pathways (Coco et al., 2013;D'Alpaos et al., 2005;Feola et al., 2005;Hughes, 2012;Zhou et al., 2014) and, hence, regulate sediment distribution within the tidal system. In addition, sediment availability determines the fate of tidal flats and salt marshes in response to increasing sea levels and wave-induced lateral erosion (Mariotti & Fagherazzi, 2013). In turn, intertidal platform elevation within the tidal frame controls the tidal prism, ultimately determining the size and extent of tidal channel networks (e.g., D'Alpaos et al., 2006;Kleinhans et al., 2015;Stefanon et al., 2012;van der Wegen et al., 2010).
Vegetation dampens tidal flow energy and strengthens the soil through the roots. Vegetation thus reduces erosional processes and enhances sediment deposition through trapping and organic soil production (Cahoon et al., 2021;D'Alpaos & Marani, 2016;Mudd et al., 2010;Temmerman et al., 2005;Toy et al., 2002). In the early stages of salt marsh formation, the flow concentration favored by pioneering vegetation patches causes erosion and channel initiation (Temmerman et al., 2007Vandenbruwaene et al., 2011;van Oyen et al., 2014). However, the reduction experienced by the tidal prism as marsh elevation increases leads to weaker in-channel currents and counteracts channel scouring . The elevation of the intertidal platform in the tidal frame, in turn, dictated by sediment availability and the rate of RSLR, controls which effect prevails over the other (Sgarabotto et al., 2021). In general, the width-to-depth ratio of salt marsh channels is smaller than that of channels cutting through tidal flats Fagherazzi et al., 2004;Lawrence et al., 2004;Marani et al., 2002;Vandenbruwaene et al., 2012;Sgarabotto et al., 2021).
RSLR adds a further degree of complexity. The growth of halophytic vegetation depends on bed elevation within the tidal frame (Morris et al., 2002;Mudd et al., 2004). Long periods of submergence are unsuitable for the survival of halophytic plants because of the poor soil aeration. High rates of RSLR and low sediment supplies threaten the development and maintenance of salt marshes. A scarce sediment input implies a weaker net deposition, leading to a smaller bed accretion and longer periods of submergence. When marsh elevation falls below a given threshold, vegetation cannot survive, and the marsh drowns (Coco et al., 2013;D'Alpaos et al., 2011;Kirwan & Temmerman, 2009;Kirwan et al., 2010;Marani et al., 2007). Marshes can thus keep pace with RSLR only if sediment supply and organic soil production are high enough to counterbalance the effects of rising sea levels Kirwan et al., 2010;Lovelock et al., 2011).
Only a few studies simulate the intertwined eco-morphodynamic evolution of salt marshes and tidal channels cutting through them (Belliard et al., 2015;D'Alpaos et al., 2007;Mariotti, 2018;Sgarabotto et al., 2021;Temmerman et al., 2007), accounting for the mutual interactions among these morphological units. These interactions are strongly influenced by the different flow resistance and erosion/deposition rates which characterize various morphological units composing an intertidal system. In this contribution, we analyzed this system from an overall perspective. We developed a mathematical model which simulates the coevolution of tidal channels and salt marshes, accounting for the higher flow resistance and the increased sediment deposition associated with the presence of vegetation, as well as for wetting-drying processes typical of shallow tidal environments. The relative importance of eco-morphodynamic feedbacks is unraveled by analyzing the evolution of an idealized tidal environment composed of narrow intertidal areas flanking a large tidal channel. Some examples of this type of tidal environments are shown in Figure 1. Typically, a series of almost evenly spaced channels, directed perpendicularly to the main channel, feed, and drain the considered intertidal area, consisting either of bare tidal flats (Figure 1a) or salt marshes (Figures 1b-1d).
Starting from an initially flat intertidal platform, we simulate the evolution of the system under the action of prescribed tidal forcing, sediment input and sea level rise. Vegetation is assumed to encroach the platform as a reference threshold elevation is exceeded. Three different types of vegetation species have been modeled separately: high marsh vegetation, low marsh vegetation, and a combination of the two mimicking the presence of coexisting species (e.g., Marani et al., 2013). Vegetation patches develop in a dynamic bed topography characterized by extending tidal channels. Vegetation feedbacks on tidal channel structure are elucidated by comparing the morphologies computed in the presence of vegetation with those obtained considering the evolution of bare tidal flats subjected to the same forcings. By synchronously simulating the intertwined evolution of salt marshes and tidal channels cutting through them, the present study shows that vegetation favors the development of channel networks with more complex structures and a higher drainage density. Under constant rates of sea level rise, a morphological equilibrium can still be attained. However, in these conditions, inner marsh portions might drown, reducing vegetation effects on channel network function.

Model Setup
The proposed eco-morphodynamic model includes four modules to simulate the interactions among the three main morphological units (tidal channels, tidal flats, and salt marshes) which generally compose a tidal basin. The different modules deal with tidal currents, sediment transport, bed evolution, and vegetation growth. Before describing in detail these modules, it is worth briefly describing how the modules are related to each other.
The structure of the numerical model is shown in Figure 2. Essentially, we assume that hydrodynamics and bed evolution are characterized by different temporal scales. In other words, the flow field is taken to adapt almost instantaneously to changes in bed elevation. These latter changes, in fact, are supposed to be relatively smooth in space and need several tidal cycles for exerting their feedback on the flow field. At each time step of a characteristic tidal cycle (e.g., a semidiurnal tide with a prescribed amplitude), the spatial distribution of flow velocity is used to solve the advection-dispersion equation for the suspended sediment concentration. At the same time, cumulative erosion/deposition rates are computed at each bed location based on the spatial distribution of bed shear stresses exerted by the flow. In the presence of vegetation, trapping of sediment and organic soil production are accounted for in the computation of deposition, while erosion is taken to be negligible. Vegetation biomass, which is computed as a function of bed elevation, modifies the friction term in the hydraulic module. At the end of any characteristic tidal cycle, the bed elevation is updated by employing the tidally averaged erosion/deposition rates multiplied by a morphological factor. The characteristic tidal cycle is then repeated.

Hydraulic Module
We consider an intertidal basin dominated by tidal forcing. As a first approximation, wind action and possible river inflows are not considered. The flow field is described through a suitably simplified version of the two-dimensional (2D) depth-averaged momentum and mass conservation equations developed by Defina (2000) to account for wetting and drying processes. At each instant t of the tidal cycle, the free surface elevation h is assumed to be given by the sum of a constant contribution ξ, coinciding with the water elevation prescribed at the seaward boundary, and a variable contribution, ζ, which expresses the local variation of the free surface with respect to the average water level (Figure 3), namely: where x is the landward directed Cartesian coordinate, with origin at the seaward boundary and y is the transverse Cartesian coordinate. Both the quantities h and ζ are referred to MSL. Owing to the relatively small flow depths typical of the investigated flows, friction is assumed to prevail over inertia and, consequently, to balance the free surface slope in the momentum equations Rinaldo et al., 1999). By evaluating the bed shear stress through the Gauckler-Strickler resistance law, it is easily demonstrated that the mass and momentum conservation equations take the form Here, U x and U y are the longitudinal and transverse components of the depth-averaged velocity, respectively; D is the effective water depth; K s is the Gauckler-Strickler resistance coefficient; ψ and ϕ are two functions (the first dimensionless) arising from the depth-averaging of the relevant equations over a representative elementary area to account for wetting and drying effects (Defina, 2000), and erfc is the complementary error function. The quantities D, ψ, and ϕ depend on the dimensionless variable where z b is the bed elevation referred to MSL and e is a characteristic scale of variations of bed elevation over the representative elementary area (i.e., a subgrid averaged roughness).
Note that the method chosen here to deal with wetting/drying processes allows one to always solve the same set of equations, without the need of removing/adding elements as the domain gets dry/wet. If the water level is much higher than the bed level, the effective water depth becomes equal to the local flow depth h − z b . Conversely, if h − z b is very small (even negative if the bed gets dry), D approaches a very small positive value that simulates a partially wet condition.
Recalling the definition of the free surface elevation (Equation 1), the simplified momentum Equations 3 and 4 can be rewritten in the form Substituting these expressions into (Equation 2), and setting = √ 2 + 2 , we obtain the following equation (van Oyen et al., 2014) that describes the spatio-temporal variation of the water level ζ(x, y, t) given the water elevation at the seaward boundary ξ.
Equation 11 is solved numerically through finite differences. The spatial derivatives are discretized through central differences, while forward differences are employed for the temporal derivative ∂ζ/∂t. The resulting system of equations is solved with the open-source computational package "Pardiso" (Schenk et al., 2001). At each time step, Picard iteration (Huang et al., 1996) is used to deal with the nonlinear terms involving U. At the first iteration, U is set equal to the magnitude of the velocity at the previous time step. After solving for ζ, at the next iteration, U is updated based on the values of U x and U y given by Equations 9 and 10, respectively. The iteration stops when the percentage difference between the values of ζ computed at previous and current iterations is less than a prescribed small value (5%).
The boundary conditions associated with the system of partial differential Equations 9-11 consist of the water level at the sea boundary, where ξ is taken equal to a cosine tidal wave with amplitude a 0 , while no-flux conditions are imposed at the other boundaries of the simulated tidal basin. The geometry of this domain is described below, in the paragraph concerning the design of numerical experiments.
At the initial instant of any simulation, a constant water level is prescribed throughout the basin at MSL, while the flow velocity is set everywhere to zero. The amplitude of the tide is then increased gradually up to the prescribed value. This relaxation procedure is employed to reduce the scouring potential in the early stages of the morphodynamic evolution and, hence, to avoid spurious numerical instabilities. The simulation of each cycle starts with the ebb phase.

Sediment Transport
Sediments are assumed to be transported mainly as suspended load, owing to their fine size. Their dynamics is described by the two-dimensional advection-dispersion equation (e.g., D'Alpaos et al., 2007) where C is the depth-averaged suspended sediment concentration, U is the depth-averaged velocity vector, k m is the horizontal mixing coefficient, and Q e and Q d are the rate of erosion and deposition, respectively. Equation 12 is solved numerically by an explicit central difference method.
During the flood phase, a given external suspended sediment concentration C sea is imposed at the seaward boundary of the tidal basin. Conversely, during the ebb phase sediments are transported out of the tidal basin. The suspended sediment concentration at the seaside boundary is thus determined by the sediment concentration coming from upstream. According to these boundary conditions, the suspended sediment is conveyed into the tidal basin and redistributed therein by the flood currents. Conversely, ebb currents usually tend to flush out the sediment. Finally, at the beginning of the simulation, the suspended sediment concentration within the basin is taken to linearly decrease from C sea at the seaward boundary to 0 g/m 3 at the landward boundary.

Bed Evolution
Local erosion and deposition drive bed evolution of the channeled and unchanneled portions of the basin. The contribution of bedload to bed evolution is assumed to be much smaller than the contribution of suspended load. Bed-elevation changes are described through the sediment balance equation (e.g., Marani et al., 2007;Toffolon & Lanzoni, 2010) where Q e and Q d are the erosion and deposition rates, respectively and R is the rate of RSLR.

of 25
The erosion flux Q e is determined on the basis of the excess of the shear stress modulus τ with respect to the critical stress τ ce . Conversely, the deposition flux in general depends on the settling velocity and on the magnitude of suspended sediment as clarified in the following.
At each location of the tidal basin the modulus of τ is computed as √ 2 + 2 , where the components τ x and τ y of the bed shear stress vector are determined by considering the local flow depth and the local free surface slope, namely The erosion rate is eventually computed as where Q e0 is a typical intensity of erosion flux, depending on the type of sediment composing the bed. Note that Equation 15 ensures a smooth transition from no-erosion to erosion conditions. This formulation, introduced by Carniello et al. (2012), prevents strong gradients of the erosion rate in a neighborhood of the critical threshold τ ce , thus ensuring a more gradual reproduction of scour processes.
The deposition rate is expressed as (D'Alpaos et al., 2007) where Q ds , Q dt , and Q do are the rates of deposition due to settling, vegetation trapping, and organic soil production, respectively. Settling deposition at any shear stress value is computed as where w s is the settling velocity of sediment particles.
On the other hand, trapping deposition and organic deposition are related to vegetation biomass, as described below.

Vegetation Growth
The local annually averaged biomass production is expressed through a fitness function , which describes the relationship between the local bed elevation z b and the biomass density  , where max is the maximum biomass density, b(z b ) is the dimensionless biomass density, f(z b ) is a fitness function and f max is its maximum value. The fitness function is species specific. Here, following Marani et al. (2013), the fitness function is described through the relation, where the elevation parameter z 0v corresponds to the optimal elevation for vegetation growth (i.e., at which the maximum biomass is attained), and the dimensional parameters, λ 1 (m −1 ) and λ 2 (m −1 ), control the range of bed elevations to which vegetation is adapted. Large values of λ are typical of specialized vegetation species, that fit within a narrow range of elevations. Conversely, small values of λ characterize species which are relatively well adapted to a broader range of marsh elevations within the tidal frame.
Vegetation patches influence the flow field by modifying the resistance to the flow. A higher biomass density leads to higher friction and, hence, a lower resistance coefficient. The overall bed friction is then assumed to be given by the sum of the local bed friction and the friction induced by the vegetation. This overall bed friction is inversely proportional to the square of a total Gauckler-Strickler resistance coefficient and is expressed as where K sb and K sv are the Gauckler-Strickler coefficients related to bed friction and vegetation, respectively.
Vegetation patches also influence bed evolution by trapping suspended sediment, through stems and leaves, and producing organic soil (the last two terms in Equation 16). Based on the approach of Mudd et al. (2004), Palmer et al. (2004), and D' Alpaos et al. (2006), these contributions to the rate of variation of bed elevation can be expressed as where Q do0 is a typical deposition rate specified empirically to account for above ground and below ground production of organic soil Mudd et al., 2004;Randerson, 1979), and ϵ v is a capture efficiency coefficient (Palmer et al., 2004). This latter coefficient can in general be related to the diameter d with α ϵ , β ϵ , and γ ϵ empirical coefficients.
Finally, the values of the diameter, the density and the average height of vegetation stems can be estimated on the basis of local vegetation biomass through the empirical relations (Mudd et al., 2004 where α n , β n , α h , β h , α d , and β d , are empirical coefficients.

Design of Numerical Experiments
In this research, we focus on the interactions between tidal channels and vegetation patches that form on a tide-dominated tidal flat of relatively small dimensions, that is, 400-m wide and 500-m long. Figures 4a-4c show the initial bed topography, assumed to be horizontal with small perturbation patches and with a base elevation of 0 m above MSL. The basin is connected with a deep and large tidal channel at the seaward boundary. A square computational grid of size 2 m is used to discretize the relevant equations.
The initial perturbation patches (in total 450) are generated on a horizontal bed and are meant to mimic the initial settlement of vegetation patches (Figures 4a and 4b). The location and magnitude of these perturbations are determined by a stochastic procedure. First, the location (x, y) of each perturbation is randomly selected. Second, a random height is imposed at each perturbation location. Then, the elevation around the selected point is averaged iteratively to generate a perturbation patch with a smoothly varying surface. This procedure is repeated a certain number of times (15 in the considered case) to expand the patch extension while reducing the patch height, eventually producing relatively isolated perturbation patches. Finally, the maximum elevation of the resulting bed topography is set equal to 0.1 m, and the height of other points is adjusted proportionally. Note that this stochastic procedure also leads to the formation of clustered patches (Figure 4b).
The water level imposed at the seaward boundary reproduces a semi-diurnal tidal cycle of period 12 hr and amplitude 0.75 m with respect to MSL (set at 0 m). A no-flux condition is instead prescribed on the other three boundaries of the basin. A sediment concentration C sea of 10 g/m 3 is imposed at seaward boundary during the flood phase of the tide, to mimic an external input of sediment. This input is in general associated with isolated recurrent events (e.g., wind resuspension, spring tides, floods in river nourished environments) but, in the long term, can be modeled imposing a constant (average) C sea as done here.
Because of the significant elevation drop that, during the ebb phase, establishes between the tidal flat and the seaward tidal channel feeding the system (see e.g., Figure 4c), extensive erosion spots may take place at the beginning of the morphodynamic evolution. In order to mitigate these localized phenomena, which may lead to numerical instabilities, the amplitude of the forcing tide is gradually increased from 0 to 0.75 m within the first 800 tidal cycles.
Typically, a control section forms nearby the seaward border of the tidal basin as the ebb level of the forcing tide becomes smaller than bed elevation (Figure 4c). This section is characterized by a Froude number F r equal to 1 and disconnects the tidal basin from the seaward channel feeding it. Given this control section, upstream (i.e., landwards) the flow is subcritical (F r < 1), while downstream (i.e., seawards) it is supercritical (F r > 1). The tide level at the seaward boundary thus no longer influences the ebb current still flowing within the tidal basin. The ebb flow is instead controlled by the condition F r = 1 that establishes in the control section. It is thus necessary to modify the seaward boundary condition, by limiting to 1 the corresponding Froude number. During the ebb phase, the effective water depth at the seaward boundary D sea is then prescribed as where the superscript i indicates the ith time step and it has been assumed, for computational convenience, that the Froude number at the boundary can be at maximum equal to 0.9. The condition (Equation 25) is retained also during the flood phase until the level of the imposed tide exceeds that of the seaward basin border. Henceforth, the flow becomes subcritical everywhere and, hence, is completely controlled by the imposed tidal wave level.
As widely adopted in long-term simulations (Coco et al., 2013;Roelvink, 2006), a morphological factor is used to speed up the morphodynamic evolution. In particular, following previous works on the morphodynamics of fluvial and estuarine environments (e.g., George et al., 2012;Morgan et al., 2020), the morphological factor is varied during a simulation, adapting its value to the overall rate of bed variation. After each simulated tidal cycle, the elevation of each computational grid point is first updated through Equation 13 and then multiplied by the morphological factor. Considering that the evolution rate is more intense during the early stages of morphodynamic evolution, we use a morphological factor of 40 for the first 100 simulation cycles (corresponding to 4,000 real tidal cycles), during which tidal channels and salt marshes are typically found to form. Later on, the erosion rate reduces significantly, and the morphological factor is increased to 200 to reproduce the longterm evolution of the system. Simulations are stopped when the system reaches an approximate equilibrium state, with almost negligible variations of the bathymetry. For the basin geometry and the set of parameters adopted in the present study, 500 simulations cycles (corresponding to 100 × 40 + 400 × 200 = 84,000 real tidal cycles, i.e., to 115 years) are always sufficient to achieve an equilibrium configuration.
Three different types of vegetation have been considered in the simulations (Figure 4d). All vegetation types start to grow as bed elevation exceeds MSL and survive up to mean high water level (MHWL). The first vegetation type (Vege #1) is better adapted to grow at high elevations, and its biomass changes with elevation similarly to that of the high vegetation species considered by Marani et al. (2013). Vege #1 biomass increases relatively slowly and reaches a maximum at 0.60 m above MSL. Then biomass decreases rapidly as bed elevation increases, eventually vanishing as MHWL is approached. The second vegetation type (Vege #2) is better adapted to grow at low elevations. This kind of vegetation mimics pioneer species, such as Spartina Alterniflora colonizing the northeastern American coast (Belliard et al., 2015;Morris et al., 2002). Its biomass production peaks at 0.15 m above MSL, and gradually vanishes toward MHWL (e.g., Marani et al., 2013). Finally, the third vegetation type (Vege #3) is taken to reproduce biomass distribution for a combination of coexisting vegetation species with optimal biomass production occurring at progressively higher bed elevations, mimicking the case of multiple vegetation species encroaching the marsh platform (Belliard et al., 2015;D'Alpaos et al., 2007;Marani et al., 2007). Just above MSL, the biomass function is similar to that of Vege #2. However, after the maximum biomass production is attained (at 0.15 m above MSL), biomass keeps constant.
It is important to note that vegetation seeds cannot settle in areas undergoing intense erosion and, hence, halophytic vegetation is significantly endangered in these areas. Therefore, in present simulations, the local biomass is set to 0 where erosion is stronger than deposition. The values of the physical and empirical parameters adopted in the various simulations are listed in Table 1. Finally, although the parameters for Vege #3 are the same as Vege #2, in the case of Vege #3, the dimensionless biomass density is kept fixed and equal to 1 (green dashed line in Figure 4d) for bed elevations higher than the optimal one (0.15 m). Using the same initial conditions, two series of simulations have been carried out by considering either a steady MSL or a MSL increasing at a constant rate. In addition, model results have been tested through a comparison with the morphological data extracted from two typical field sites located in the Venice Lagoon (Italy) and in Warbah Island (Kuwait).
During the morphodynamic evolution of the basin, tidal channels cut through both tidal flats and salt marshes. The corresponding bankfull channelized areas are distinguished using the channel detection method proposed by Geng et al. (2018). This method identifies channelized areas on the basis of the local variations of bed surface elevation and has been demonstrated to be fairly robust in delineating the structure of creek networks in intertidal zones with complex topographies and rapidly varying flow field.

Results
The initiation and development of tidal channels and the accretion of the intertidal platform are simulated for bare soil conditions and for each of the three vegetation types previously introduced. In the following sections, we summarize the various morphodynamic processes emerging from the simulations.

Coupled Evolution of Salt Marshes and Tidal Channels
At the beginning of each simulation, the tidal amplitude is gradually increasing and has not yet reached its final value. The resulting flow field is relatively weak and hence deposition and erosion rates are quite small. As the tidal amplitude approaches its final value (after 800 tidal cycles, i.e., about 1.1 yr), erosional processes amongst bed perturbation patches gradually increase. Tidal channels start to grow from the seaward boundary through headward erosion after about 2.5 years and a still immature channel network is visible after about 4 years (Figures 5a, 5e, 5i, and 5m). The structure of this incipient network consists of three main channels with similar spacing and, hence, similar drainage areas. This structure is quite similar for all the simulations, independently of the presence and type of vegetation. Indeed, after 4 years the elevation of the intertidal platform is still too close to MSL and, when vegetation is considered, the biomass is too small to affect significantly channel network morphology and structure.
As the sediment enters the tidal basin through the seaward boundary, the elevation of the intertidal platform progressively increases. The suspended sediment concentration decreases landward, owing to deposition induced by settling and by a progressive decrease of advective transport as prescribed by the advection-diffusion Equation 12. Before the formation of channels, relatively high velocities (order of 0.3 m/s during the ebb phase) occur on the tidal platform near the seaward elevation drop, prompting erosion rates higher than deposition ones. Conversely, proceeding landward, velocities become weaker (order of 0.2 m/s during the flood phase) and, hence, deposition overcomes erosion.
After their initial growth, channels keep deepening and extending landward, leading to more complex channel networks. Perturbation patches prompt channels to curve and to branch diverting or splitting their flows, respectively ( Figure 5). This general behavior is enhanced by vegetation. After 12 years, channel branching and lengthening occur more frequently in the case of Vege #2 and Vege #3, whose biomass can attain high values also during the early stages of the morphodynamic evolution (Figures 6f-6h, and 6j-6l).
Conversely, in the case of Vege #1, the optimal bed elevation for biomass growth is close to MHWL (Figure 6d). At this marsh elevation, the channel network has already reached a relatively defined configuration. Consequently, the differences with respect to the channel network obtained without vegetation are less evident (Figures 5b-5d, and 5f-5h). Instead of extending from the seaward boundary, some isolated channel portions form directly in the middle of the tidal basin, especially in cases of Vege #2 and #3 (Figures 5j, k and 5n, 5o). These isolated channels are located between vegetation patches (Figures 6f and 6j), and their inception is driven by the strong flow concentration and the corresponding erosion taking place there. Note that in all four cases, a channel generates and extends landward along the lower boundary (y = 500 m). This is due to the structure of the initial randomly distributed bed perturbation that enhances flow concentration near the lower boundary, triggering the formation of the straight channel.
During the bed evolution, platform elevation controls vegetation encroachment and, through the adopted biomass parametrization, its effects on morphology. In the case of Vege #1, as the mean bed elevation of the tidal basin gradually increases, vegetation first starts to encroach seaward areas and then extends landward (Figures 6b-6d). This trend is qualitatively similar to that observed in Figure 1b. In the cases of Vege #2 and Vege #3, biomass production is quite high also during the early stages of evolution, when bed elevations are still close to MSL. This allows salt marshes to initially form in the middle and landward areas of the tidal basin, where deposition rates exceed erosion rates. As the intertidal platform emerges on average and exceeds the optimal elevation for biomass, Vege #2, progressively decays (Figures 6f-6h). Conversely, Vege #3 gradually spreads throughout the tidal basin (Figures 6j-6l) enhancing deposition rates. At the end of the evolution, the mean elevation of the salt marshes is thus higher compared to the other cases.

Overall Vegetation Effects
The overall effects of vegetation on the various morphological units, whether channels, marshes or tidal flats, are here evaluated. To this goal, we tracked the temporal evolution of the total vegetation biomass,  , the total channel length, L c , the volume of the channel network, V c , the mean length of unchanneled flow paths, ℓ, the drainage efficiency ℓ H /ℓ, and the cumulative amount of erosion  . The total channel length, L c , is computed by summing the axis length of all the channels cutting through the tidal basin. The total channel volume is defined as the volume between the channel bed and the elevation of the intertidal platform at the channel borders (Geng et al., 2018). The unchanneled flow path length, ℓ, is computed as the mean distance from a point on the intertidal platform to the nearest channel (Marani et al., 2003). The drainage efficiency is given by the ratio ℓ H /ℓ, with ℓ H the Hortonian path length given by the ratio of tidal basin area to total channel length (Horton, 1945). For a given Hortonian length, ℓ H , high values of ℓ H /ℓ correspond to small values of ℓ, indicating that the spatial arrangement of the tidal channels efficiently reduces the mean overmarsh path length (Marani et al., 2003). Finally, the cumulative erosion  is computed as the total amount of sediment eroded during the evolution.
All the above defined global variables are plotted as a function of time in Figure 7. Both the channel length and the volume of the channel network grow rapidly at the beginning of the evolution (Figures 7a and 7b). After about 14 years, as the channels have almost extended throughout the whole tidal basin, the rate of channel lengthening slows down ( Figure 7a). As the channels grow landwards, the mean unchanneled length decreases rapidly, eventually tending to a constant value after the total channel length has approached its maximum value (Figure 7c).
The development of vegetation not only promotes channel lengthening (Figure 7a) but also reduces the overmarsh unchanneled length ( Figure 7c). As compared to bare soil conditions, the drainage efficiency of the tidal network increases for vegetation #2 and #3, while it decreases for vegetation #1 (Figure 7d). These findings can be explained by considering the spatial distributions of bed elevation and the channel morphology shown in Figure 5. It emerges that in the case of Vege #1, the final distribution of biomass ( Figure 6d) favors channel formation in the seaward and middle portions of the tidal basin, leading to more channel branches than in the bare soil case (Figure 5h). Conversely, the most landward zones are characterized by a lower degree of channelization, thus leading to a reduction of the overall drainage efficiency of the tidal basin. Figure 6. The spatial distribution of vegetation biomass (g/m 2 ) is plotted at different evolution stages (4, 12, 45, and 115 years) for the three vegetated cases. Vege #1, Panels (a-d); Vege #2, panels (e-h); Vege #3, panels (i-l).
In the cases of Vege #2 and Vege #3, biomass initially reaches quite high values throughout the entire tidal basin (e.g., Figures 6f and 6j). The tidal channels develop all over the basin and develop many extra branches. Consequently, the drainage efficiency raises (Figure 7d). In particular, the most efficient drainage system is obtained in the case of Vege #3 (Figure 5p) which, at the end of the simulation, exhibits the maximum total biomass as well (Figure 6l). The maximum biomass is approximately four times larger than in the case of Vege #2 (Figure 7f). Indeed, in this latter case, vegetation biomass progressively decreases as the bed elevation further accretes after exceeding the optimal level for biomass production. Note that the jumps in drainage efficiency observed for Vege #3 at about 65 and 107 years are associated with the simultaneous headward growth of various channels in the innermost portions of the tidal basin.
It is also worthwhile to note that, during the early stages of evolution (27 years), cumulative erosion in the cases of Vege #2 and Vege #3 increases much faster than in the other two cases (Figure 7e). Afterward, differences in cumulative erosion become increasingly smaller and, at the end of the various simulations, tend to vanish. This finding indicates that, even though the final equilibrium configurations of the tidal basin can be characterized by different morphologies (Figure 5), they are obtained with almost the same total cumulative erosion.  Figure 8a). Later on, the mean elevation of unchanneled areas becomes higher when vegetation patches spread throughout the entire basin and the biomass is such to yield a significant production of organic soil (Vege #2 and #3, Figures 8d and 8g). Conversely, for slowly growing vegetation patches (Vege #1) the differences with respect to the bare soil case are fairly small. Moreover, in the absence of vegetation, the landward portion of the platform is characterized by slightly higher elevations as compared to the case of Vege #1. This occurs because vegetation patches in the seaward and middle portion of the basin trap the suspended sediment coming from the sea, and thus reduce the amount of sediment delivered landwards leading to lower deposition rates. This overall behavior is consistent with field measurements carried out by Gong et al. (2017) along a cross-shore profile of the Jiangsu coast (China). The temporal variations of bed elevation recorded along this profile during about one year indicate that cumulated bed accretion attains a maximum in the transitional area between bare tidal flats and densely vegetated salt marshes. Proceeding landwards, sediment trapping by vegetation leads to a decrease in cumulated bed accretion, and a slight bed erosion is observed at the gauging station located in the innermost vegetated areas . In the present simulated tidal basin, the final equilibrium configuration of the intertidal platform exhibits higher Figure 8. The longitudinal distribution of the mean bed elevation (m above MSL) in channelized and unchannelized areas (panels a, d, g), along-channel distribution of cross-sectional width (panels b,e,h), and cross-sectional maximum depth (panels c, f, i) are plotted at different times of the morphodynamic evolution: 12 years, panels (a-c); 45 years, panels (d-f); 115 years, panels (g-i). The cross-sectional parameters are measured from a typical channel located in the middle of the tidal basis. Black dots refer to bare soil while red, blue, and green dots denote Vege #1, Vege #2, and Vege #3, respectively. elevations near the sea boundary, it decreases progressively landwards for the first 120 m and then it keeps almost constant or slowly decreasing up to the innermost areas.
The overall morphological variations of channel cross sections can be characterized through the bankfull cross-sectional width and the maximum cross-sectional depth. The along-channel distribution of these quantities is plotted in Figures 8b, 8e and 8h and c,f,i with reference to a representative main channel located in the middle of the tidal basin (y = 200-300 m). During its development, the channel progressively extends from the seaward border to inner areas of the tidal basin, developing mild bends and minor branches.
A slight narrowing invariably characterizes the cross sections close to the seaward boundary, where the adjacent unchanneled areas are dominated by deposition. The cross-sectional width then increases as the channel extends in the central region of the tidal basin and further decreases toward the landward channel head (Figures 8b, 8e and 8h).
In the considered tidal basin, higher elevations of intertidal areas boosted by vegetation growth enhance channel widening. This effect is particularly evident in simulations carried out with Vege #2 and Vege #3 (blue and green dots in Figures 8e and 8h). In these two simulations, as compared to Vege #1 and bare soil cases, sediment trapping and organic soil production lead to higher bed elevations of unchanneled areas in the central and landward portions of the basin (blue and green continuous lines in Figures 8d and 8g). On the contrary, near the seaward border, similar channel widths are attained in all four cases. This is due to the intense settling deposition that takes place in this area, which definitely prevails over the contribution of vegetation to sediment trapping and organic accretion.
Channel depth increases over time approaching an almost equilibrium condition as the rate of change in bed elevation tends to vanish. During the first 12 years, the channel deepens quite rapidly while extending landwards (Figure 8c). After ∼45 years, the channel depth is quite similar for all simulations (Figure 8f). The presence of vegetation eventually (after 115 years) leads to shallower cross sections in both the seaward and central areas of the basin (Figure 8i). This effect is particularly evident for Vege #2 and Vege #3. The higher mean elevation of the adjacent intertidal platform characterizing these simulations, as compared to the bare soil case, implies that a smaller volume of water needs to be accommodated in channels during the ebb tide. Given also the larger cross-sectional width in the case of Vege #2 and #3 (Figure 8h, blue and green dots), this reduction in the ebb flow discharge explains why the seaward and central reaches of the channel get shallower (Figure 8i). The larger channel widths and shallower depths in the vegetated cases lead to width-to-depth ratio higher than in the bare soil situation, differently from field observations Lawrence et al., 2004;Marani et al., 2002). This discrepancy is largely due to the lack of a suitable subgrid model of bank erosion. In addition, the method used to detect channelized area tends to overestimate the channel width in the presence of vegetation as compared to the bare soil case. Indeed, bare soil channels exhibit a more delineated channel border with respect to vegetated soil. For vegetated cases, the stronger deposition near the channel banks leads to a smoother transition of the cross-sectional profile from the steep channel bank to the nearly horizontal intertidal platform adjacent to it.

Effects of Sea Level Rise
In coastal areas, intertidal zones are usually influenced by increasing MSL. Mathematical models provide a fundamental tool to evaluate the possible long-term consequences of changes in MSL on tidal eco-morphodynamics. Different from well-developed salt marshes which have been widely studied (e.g., Lovelock et al., 2011;Rodríguez et al., 2017), in this study, we focus on channel-marsh interactions from the initiation phase to the final equilibrium configuration. In order to unravel the complex feedbacks between morphodynamics, vegetation dynamics, and environmental forcing, we have employed constant rates of RSLR. In the case of the schematic tidal basin considered here, keeping fixed all the parameters listed in Table 1, no major differences are observed for relatively low rates of RSLR (2 mm/yr). Conversely, under a fairly large rate of RSLR (8 mm/yr), platform elevations fall below MSL, thus creating unsuitable conditions for vegetation growth. In the following, we discuss the eco-morphodynamic changes experienced by the investigated tidal basin when subject to a rate of RSLR of 4 mm/yr, close to the value of 3.5 mm/yr usually adopted for the Venice Lagoon in long-term simulations (e.g., Marani et al., 2010).
The bed evolution of the four scenarios under RSLR is shown in Figure 9. In order to synchronously denote salt-marsh development, vegetation patches with biomass density larger than 100 g/m 2 are denoted by dashed lines. Similar to simulations carried out in the absence of RSLR, tidal channels start to grow from the seaward boundary of the tidal basin and then extend landward through headward erosion. After channel formation, the mean basin elevation gradually increases due to a positive net deposition rate. Vege #2 and Vege #3 patches start to grow throughout the entire basin from the very beginning of the evolution (areas denoted by black dashed lines in Figures 9i-9l and 9m-9p), while Vege #1 start to develop at higher elevations and, hence, begins to encroach the intertidal platform later (Figures 9e-9h). However, RSLR dramatically restricts vegetation growth in the inner areas of the tidal basin as compared with constant MSL simulations. Indeed, because of the progressive landward decay of the transported sediment, these areas Figure 9. The spatial distribution of bed elevation (m), referred to mean sea level (MSL), is plotted at different evolution stages (4, 12, 45, and 115 years) for a sea level rise of 4 mm/yr. The four simulated cases are the same of Figure 5: bare soil, panels (a-d); Vege #1, panels (e-h); Vege #2, panels (i-l); Vege #3, panels (m-p). Black continuous lines denote the edge of tidal channels. Black dashed lines delimit vegetated areas with biomass density larger than 100 g/m 2 .
hardly keep pace with the rate of RSLR and become less prone to vegetation growth. Only some vegetation patches with significant biomass survive because of the locally enhanced bed accretion induced by sediment trapping and organic soil production. Outside of these patches, bed accretion is definitely lower and bed accretion cannot keep pace with the rate of RSLR, which threatens marsh survival. Highly vegetated patches thus hardly extend and become increasingly isolated. This uneven spatial distribution of vegetation exacerbates the differences in topography across the tidal basin. Eventually, vegetation biomass inevitably concentrates near the seaward border, where higher bed elevations are attained. The basin thus develops a remarkable landward-decreasing bed slope characterizing the basin in these scenarios (Figures 9d, 9h, 9l, and 9p).
The overall effects of RSLR on tidal-channel and marsh morphology are shown in Figure 10 in terms of total channel length and total biomass. In the unvegetated case, the increasing mean sea level has little influence on the development of the tidal channel networks, at least for the present schematic basin. As compared with the case with no RSLR (Figure 7a), the total channel length measured at the end of the simulation increases slightly, by about 2.1%.
On the contrary, RSLR has a strong impact on total vegetation growth. Indeed, owing to the overall lower relative bed elevation, vegetation biomass is remarkably smaller than in simulations carried out with a constant MSL (Figure 7f). As a result, vegetation affects the development of tidal channels to a less extent. The total biomass of Vege #1 is an order of magnitude smaller than that of Vege #2 and Vege #3 (red dotted line in Figure 10). Vegetation growth concentrates near the seaward border of the basin where, as shown by the dashed lines in Figure 9h, higher bed elevations ultimately occur. The total channel length is about 4% shorter than that observed in the case of a constant MSL (see the red continuous line in Figure 10 as compared to Figure 7a).
Also for the other two vegetation types, RSLR causes a remarkable decrease of total biomass, which at the end of the simulations turns out to be approximately six times (Vege #2) and four times (Vege #3) smaller than in the case of a constant MSL (see the blue and green dotted lines in Figure 10 as compared to Figure 7f). Total biomass reduction is particularly severe in the period between 14 and 34 years when vegetation patches, initially grown in the inner portions of the tidal basin, progressively decrease in size owing to the reduced accretion rate of the intertidal platform. In both cases, similarly to Vege #1, biomass eventually concentrates near the seaward basin border (the dashed lines in Figures 9l and 9p). As compared to the bare soil situation, the influence of Vege #2 and Vege #3 on the tidal channel morphology is fairly similar and prompts a slight increase (about 43%) in the total channel length (blue and green dotted lines in Figure 10).

Comparison With Field Data
Multiple environmental factors influence the development of tidal channel networks, leading to a wide range of structures and patterns. In the absence of detailed information on the boundary conditions and on the past changes in landforming processes occurred within an intertidal environment, it is almost impossible to reproduce numerically the actual morphology of a tidal channel system. A general comparison between the features of real and synthetic channel patterns may however be based on some relevant statistics (e.g., the mean channel width and length, and the mean channel spacing).
As previously stated, in the present contribution we focus on small tidal environments connecting with large channels as those shown in Figure 1, concerning the south edge of Warbah Island (Shatt el-Arab Estuary, Kuwait), the Great Ouse River (The Wash, UK), the Petaluma River (California, USA), and the Dell'Ancora Channel (Venice Lagoon, IT). Note that, while Warbah Island is characterized by channels cutting through bare tidal flats, the other three cases instead refer to salt-marsh channels.
In particular, we applied the model to mimic the features of tidal channels forming in the rectangular areas depicted with red lines in Figures 1a and 1b. The tidal basin considered in Warbah Island (hereafter denoted with WI) has a length of 600 m and a width of 500 m. In Kuwait Bay, the tides are semidiurnal and the tidal range changes from 4.2 m (during spring tides) to 0.5 m (during neap tides) (Baby, 2011). Considering the cosine tide wave used in the simulation, the tidal range is set to 2 m. Given the lack of specific information, the seaward sediment concentration is tentatively set to 10 g/m 3 , corresponding to the lower SSC values observed in the northwest Persian Gulf (Al-Ghadban, 2004;Al-Yamani et al., 2004). In the case of the salt marsh flanking the Dell'Ancora channel in the Venice lagoon (hereafter denoted as VE), the tidal basin length is 80 m and the width is 50 m. Similar to other studies carried out for the Venice Lagoon (D'Alpaos & Marani, 2016), the tidal forcing is taken to be sinusoidal with amplitude of 0.5 m. The input sediment concentration is set to 3 g/m 3 , and kept constant in time. This value has been taken lower than the sediment concentration measured in the neighboring area (about 8.6 g/m 3 , Carniello et al., 2012;Venier et al., 2014) to account for the intermittency of the wind-events leading to the observed concentration value. Two vegetation scenarios are considered. Typically, salt marshes in the Venice Lagoon are colonized by multiple vegetation species Silvestri et al., 2005). Vegetation biomass then attains high values even when bed elevation reaches relatively high values. Here, for sensitivity analysis purposes, two different biomass functions (Vege #VE1 and Vege #VE2) have been considered to mimic a multiple species vegetation ( Figure 11d). In both cases, the biomass is kept fixed after its maximum is attained for a prescribed elevation (0.4 m above MSL for Vege #VE1, and 0.1 m above MSL for Vege #VE2). Simulations lasted 115 years, such that the channels and the surrounding intertidal platform have approximately reached a stable configuration characterized by a vanishing small rate of morphological changes. Figures 11a-11c show the bed topographies at the end of the simulations. Almost parallel channels grow within the simulated basins, as observed in the field (Figures 11e and 11f). However, compared with the actual channels the numerically generated channels have fewer branches. This finding is possibly related to the grid size (2 m for simulation WI and 0.5 m for simulations VE) and the initial bed perturbations used in the simulations, which limit the formation of smaller branches.
In the case of the VE marsh, the mean width, the length, and the spacing (measured at the seaward border) of the various simulated channels are definitely similar to the actual ones ( Figure 12). Channel features strongly depend on the size of the tidal basin. Indeed, the Venice channels investigated here are definitely narrower than synthetic channels recently reproduced in Sgarabotto et al. (2021), where mean channel width reaches 20 m in a basin 100 times larger than that here considered. In the case of WI, simulated channels are generally wider and have a smaller spacing as compared to channels observed in the field. In general, the spacing between channels observed in the field increases with channel width and length (gray dots in Figure 12b). This is not surprising, as longer and wider channels feed and drain wider intertidal areas, and hence, are characterized by larger spacing. The simulations carried out for the VE marsh (red and blue crosses in Figure 12) correctly reproduce this trend. On the contrary, in the WI simulation, many channels with different lengths have similar spacing intervals, which vary in a narrow range (black cross in Figure 12b). This finding might be due to the uncertain information about model parameters in the case of WI. Note that in the simulations the channel length is likely related to the imposed basin length, which exerts some control on the volume of water drained during every ebb tide.
Clearly, various uncertainties affect the parameters used in the simulations, especially for WI (e.g., the externally imposed concentration, the critical shear stress for erosion, initial bed perturbations). These uncertainties might in part explain why the simulated WI channels exhibit a smaller range of width values as compared to reality. Certainly, the simulated environmental conditions are simpler and more homogeneous than in the field. Besides, near bank flow, and hence bank erosion, are poorly represented in the model. The initial bed topography can also affect the channel network morphology. The initial randomly generated disturbances superposed to an otherwise horizontal bed may influence the final channel network geometry, leading to less spaced channels than observed ones in the case of WI.

Model Limitations
The present modeling framework is based on some relevant assumptions. First, friction is assumed to dominate over inertia (Rinaldo et al., 1999) in the momentum equations. Nevertheless, inertia can play an important role within the channels as well as near vegetation patches, leading to longer wakes behind them (van Oyen et al., 2014). Second, erosion at the channel banks and at the seaward border of the tidal basin is described approximately as a continuous process. No subgrid parameterization is used to account for the actual shape of the bank/border and of localized and intermittent bank collapse events. These approximations surely influence channel width computations. Moreover, condition (Equation 25) imposed at the seaward boundary to limit the Froude number during late ebb and early flood may have some influence on the morphology eventually attained by the seaward border of the tidal basin. Wave-induced erosion, not accounted for in the model, may also matter even in the presence of a vegetated platform (e.g., Leonardi & Fagherazzi, 2014;Leonardi et al., 2016;Marani et al., 2011;Mariotti & Fagherazzi, 2013). In our simulations, we do not account for the effects of fluvial water discharge. The corresponding sediment input would increase sediment availability, thus favoring vertical accretion and horizontal progradation of salt marshes and, hence, their resilience to sea-level rise (Alizad et al., 2018;D'Alpaos et al., 2011;Sgarabotto et al., 2021). In addition, the present model investigates channel and salt-marsh initiation under constant rates of RSLR. Evaluating the response of these systems and their resilience to the accelerating rates of RSLR suggested by recent projections (Church et al., 2013;Oppenheimer et al., 2019) is undoubtedly an important issue (e.g., Alizad, Hagen, D'Alpaos et al., 2011;Kirwan & Temmerman, 2009;Rodríguez et al., 2017) which can be efficiently tackled through the present modeling approach. It is also worthwhile to note that the use of increasing rates of RSLR might prevent the system to reach equilibrium conditions. Finally, the offline technique for the bed update may affect the synthetic morphologies.
In the present study, we focus on the effects of different types of vegetation which singly encroach the marsh. Therefore, the present model does not account for the competition among different vegetation species and the possible change in the dominant species associated with the effects of RSLR Marani et al., 2013). Salinity dynamics and its effect on vegetation are not included in the model as well, and deserve attention in the future. In principle, the present model can be modified to simulate the eco-morphodynamic evolution of mangrove systems, which occupy a physiographic niche similar to salt marshes (Cahoon et al., 2021). Clearly, the parametrizations concerning flow resistance, sediment trapping, and organic soil production need to be suitably adapted, given the different physiology and morphology of mangroves as compared to marsh vegetation (Phillips et al., 2017;Xie et al., 2020). Overall, despite the simplifications introduced with respect to 2-D morphodynamic models (Boelens et al., 2018;Coco et al., 2013;Hibma et al., 2003;van der Wegen & Roelvink, 2008), the present model is deemed to reproduce correctly the eco-morphodynamic evolution of the various morphological units composing small marine-dominated tidal basins, producing tidal morphologies which reasonably resemble those observed in the field.

Conclusions
This study focuses on the eco-morphodynamic coevolution of tidal channels and salt marshes in a schematic basin, mimicking tidal environments flanking large tidal channels or tidal rivers. The coevolution is simulated through an eco-morphodynamic model which accounts for wetting-and-drying processes and vegetation-induced roughness. Different types of vegetation have been considered, each characterized by a specific biomass distribution depending on bed elevation. Simulations have been carried out by letting the system evolve starting from a horizontal, slightly and randomly perturbed bed, under the influence of a forcing tide and a sediment concentration imposed at the seaward boundary. The main results of our analysis can be summarized as follows.
Model results successfully mimic the morphologies observed in the field for both bare soil and vegetated conditions. In general, widely spread vegetation patches are found to definitely affect the structure of the tidal channel network. The increased friction produced by vegetation feeds back on the flow field and, consequently, on the erosion and deposition patterns that ultimately determine the structure of the tidal channel network.
Vegetation with an optimal elevation for biomass production close to MSL is found to strongly promote the development of tidal channels as compared with the bare soil case. Specifically, channels extend more rapidly landward, through headward erosion. Their cross sections are usually wider and shallower. In addition, more channel branches grow, forming more complex network structures with a higher drainage efficiency with respect to the bare soil case, as well as to the case of vegetation with a higher optimal elevation for biomass production.
In general, vegetation starts to colonize the intertidal platform from the seaward border, where the externally supplied sediment concentration ensures a more rapid bed accretion. Vegetation then extends landwards. The presence of pioneer vegetation (mimicked by introducing randomly distributed perturbations of the initial bed topography) enhances the deposition on unchanneled areas, leading to faster growth of salt marshes, especially for species with optimal biomass production occurring closer to MSL. Vegetation trapping of suspended sediment in the seaward and middle portions of the tidal basin invariably reduces the amount of sediment delivered to landward areas, weakening the sedimentation there.
In the presence of RSLR, owing to the lower relative bed elevation above MSL, vegetation growth is limited and, consequently, its control on channel form and function is reduced. Given a sufficient seaward sediment supply, the deposition rate at the seaward border of the basin can keep pace with the rate of RSLR. Conversely, the inner basin becomes incrementally submerged and a landward bed slope forms. Vegetation patches with high biomass become increasingly isolated. Indeed, the bare tidal platform adjacent to these patches experiences a deposition rate much lower than that needed to counteract the effect of RSLR.
The simulated channel networks exhibit a reasonable similarity with the parallel channel patterns observed in tidal areas adjacent to larger tidal channels. Nevertheless, the simulated cross sections are somewhat wider than those observed in the field, especially the smaller channel branches. This finding is strictly related to the relatively large grid size used in the simulations, which prevents reproducing correctly the smaller creek geometry. Clearly, the lack of a sub-grid parametrization of bank erosion and bank collapse events can also explain the differences between simulated and real channel geometries, calling for the inclusion of bank processes in future model developments.

Data Availability Statement
The data used in the current study are available at https://doi.org/10.5281/zenodo.5603005.