Effects of Vegetation, Sediment Supply and Sea Level Rise on the Morphodynamic Evolution of Tidal Channels

Tidal channels play a leading role in the morphodynamic evolution of tidal environments by driving the exchange of water, sediments, and nutrients within these systems. Despite their fundamental function on landscape construction, current knowledge on their plano‐altimetric evolution in response to environmental forcing is still limited. We developed a mathematical model to address the interplay of biogeomorphic processes on the plano‐altimetric equilibrium of a straight tidal channel flanked by adjacent intertidal platforms. The model couples a one‐dimensional hydrodynamic model with a bed evolution model. The former describes the flow field within the channel and the storage contribution of lateral shoals. The latter accounts for the main physical processes responsible for shaping the tidal landscape, namely erosion, deposition, and sea level rise, mediated by halophytic vegetation growth. Model results reproduce a variety of channel morphologies, ranging from those typical of tidal‐flat channels, with large width‐to‐depth ratios, to those characterizing salt‐marsh creeks, which have much narrower cross sections and small width‐to‐depth ratios. Our results show that vegetation encroachment on the marsh surface produces two competing effects. First, vegetation increases flow concentration within the channel, owing to the increased flow resistance on the vegetated platform thus favoring channel incision. Second, vegetation promotes a reduction in the tidal prism, owing to higher accretion rates of the salt‐marsh surface thus leading to channel infilling. Which effect prevails over the other depends on the elevation of the intertidal platform in the tidal frame. The latter is in turn dictated by sediment availability, the rate of relative sea level rise, and vegetation biogeomorphic effects.

cal features of a tidal basin, and their intertwined dynamics, are jointly tackled. The study deals with this relevant issue.
Given their fundamental control on transport processes, tidal channel morphodynamics has been addressed by many authors starting from the landmark studies by Friedrichs (1995), Schuttelaars and de Swart (1996), and Lanzoni and Seminara (2002). Friedrichs (1995) related the equilibrium cross-sectional geometry of sheltered tidal channels to the shear stress just necessary to maintain a null gradient in net along-channel sediment transport. A lower bound on this stress is provided by the critical shear stress for initiation of sediment motion. The application of this approach allowed Friedrichs (1995) to explain the classical tidal prism-channel area relationship (O'Brien, 1931;D'Alpaos et al., 2010;Escoffier, 1940;Jarrett, 1976;van der Wegen et al., 2010).
The morphodynamic pattern formation in constant-width short embayments (i.e., with length much shorter than the tidal wavelength) was instead investigated by Schuttelaars and de Swart (1996). These authors developed a model accounting for net sediment transport due to flow asymmetry, settling lag, diffusion, bed slope, and linearized friction. Model results showed that short tidal embayments attain a unique cross-sectionally averaged equilibrium state defined by a vanishing net sediment transport within the embayment. The synthetic profiles compared favorably with those observed in various embayments of the western Wadden Sea. Schuttelaars and de Swart (2000) extended the above-described modeling framework by considering embayments with arbitrary length and by applying continuation methods. Starting from the short-embayment equilibrium, they progressively increased the embayment length, locally imposing a vanishing tidally averaged sediment flux. Unique morphodynamic equilibrium configurations were found to exist only in the absence of external overtides and for embayment lengths shorter than the frictional length scale of the tidal wave. Multiple equilibrium configurations were instead obtained in the presence of significant external overtides and for a certain range of embayment lengths. The modeled equilibrium profiles were in good agreement with those observed for various short embayments of the Wadden Sea and, in the long-embayment case, for the Western Scheldt. Lanzoni and Seminara (2002) relaxed the assumption of a constant-width embayment by considering the morphodynamic equilibrium of funnel-shaped tidal channels. They considered rectangular cross sections with landward exponentially decreasing widths, fixed channel banks, and cohesionless bottom sediments. The role of channel convergence on the bed profile was addressed numerically solving the full set of one-dimensional de Saint Venant and Exner equations. A sinusoidal tide was imposed at the seaward boundary and the effects of water storage over the intertidal platforms adjacent to the channel were neglected. The equilibrium profiles were characterized by a concavity that increased with the degree of channel convergence. The associated tidally averaged sediment flux tended to vanish everywhere and a nearly constant value of the maximum flood/ebb speed was reached along the channel. The synthetic bottom profiles nicely agreed with those of various channels in the Venice Lagoon and with the experimental results by Tambroni et al. (2005).
The above results were substantially confirmed by the one-dimensional simulations carried out by Todeschini et al. (2008). These authors explored the role of the forcing tidal amplitude, bed friction, and degree of convergence on the long-term morphodynamics of tidal channels. Model results highlight the crucial role exerted by the degree of width convergence on bottom profiles, with strongly convergent channels being shorter than weakly convergent ones. This latter finding was later on confirmed by the perturbative approaches put forward by Seminara et al. (2010) and Toffolon and Lanzoni (2010). Both of these studies emphasized that, given the sediment characteristics and the tidal forcing at the inlet, channel convergence not only leads to a concave equilibrium bed profile but also tends to shorten the equilibrium channel length. In addition, the storage action exerted by tidal flats and/or salt marshes flanking the channel (an effect neglected by previous studies) was found to produce a significant shortening of the equilibrium length.
All the above recalled modeling frameworks Schuttelaars & de Swart, 1996, 2000Todeschini et al., 2008;) imposed a rectangular shape of channel cross sections and a priori longitudinal variations in channel width, or absence thereof, neglecting the mutually related adjustments between channel width and depth. Mathematical frameworks describing the three-dimensional SGARABOTTO ET AL.
10.1029/2020WR028577 equilibrium morphology of a tidal channel and of the adjacent intertidal platforms without imposing constraints on cross-section shape or size, and longitudinal variations in channel width or depth are rare (Canestrelli et al., 2007;van der Wegen et al., , 2010. These studies employed two-dimensional, process-based numerical models coupling flow and bed level changes, with large computational efforts. Exploring the influence of environmental conditions (e.g., sediment supply, sea level rise) on landscape morphodynamic evolution by using the latter modeling frameworks, would further increase the computational burden involved, thus calling for simplified modeling approaches of the type proposed herein.
The present work aims to develop a mathematical model describing the intertwined biomorphodynamic evolution and plano-altimetric equilibrium configuration of a straight tidal channel and of the adjacent intertidal platforms (see Figure 1). The considered scheme consists of a rectangular tidal basin, cut through by a tidal channel flanked by intertidal areas. The role of environmental forcing (e.g., sea level rise and sediment supply) on the evolution of the tidal landscape is accounted for, together with the crucial effects of vegetation on the dynamics of the intertidal platform and of the channel which dissects it Temmerman et al., 2007).
The starting point of our analysis is the purely erosional morphodynamic model of Lanzoni and D'Alpaos (2015) (hereafter denoted as LD15). The model, which is valid for short channels, has recently been revisited by Xu et al. (2019) to account for inorganic sediment deposition driven by an external sediment supply. In both studies, the rectangular tidal basin is discretized by a sequence of transects describing the topography of the tidal channel crosssections and the adjacent intertidal areas. Proceeding seaward, each transect feeds or drains a larger basin surface. Under the assumption of a short tidal basin, flood and ebb tidal flows are described by a quasistatic hydrodynamic approach, whereby the free surface keeps horizontal while its elevation changes through the tidal cycle. The various transects are thus linked together in terms of discharge, but not in terms of momentum. Given the discharge pertaining to a specific transect, the transverse distribution of bed shear stresses is computed through the procedure adopted for a single cross section by Furbish (2001) andD'Alpaos et al. (2006).
In the present study, we relax some of the assumptions embedded in the model by Lanzoni and D'Alpaos (2015). First, the various transects are linked together in terms of momentum and not just in terms of SGARABOTTO ET AL.  mass conservation. This allows us to overtake, from a hydrodynamic point of view, the short-channel limitation. Moreover, wetting and drying processes are now accounted for the mass conservation and momentum equations. Second, we account for the effects of an external sediment supply within the tidal basin, thus modeling the possible aggradation of the intertidal platform flanking the channel. Third, the effects of vegetation growth on the marsh platform exceeding a threshold elevation are accounted for.
The study is organized as follows. Section 2 describes the biomorphodynamic model and its various components. Section 3 presents the results of the simulations carried out to fully explore the parameter space under different scenarios of sea level rise and sediment supply, with particular emphasis on the mutual role exerted by physical and biological processes. Section 4 discusses the relevant features of the synthetic morphologies, with specific attention to the role of vegetation, compares model results with field observations, and highlights model limitations. Finally, Section 5 draws a set of conclusions.

Methods
The proposed model describes the biomorphodynamic evolution of a schematic rectangular basin of width B T and length L (  The model is composed of three modules (Figure 3). The hydrodynamic module computes the one-dimensional (1D) hydrodynamics within the tidal basin. The shear stress distribution module uses the friction slope computed for each transect to determine the bed shear-stress distribution within each transect. In the morphodynamic module, such distribution is used to compute accretion/erosion rates and determine the bed temporal evolution within the tidal basin.

Hydrodynamics
The hydrodynamics within the tidal basin is described using a cross-sectionally averaged model, whose continuity and momentum equations read (Cunge et al., 1980;Dronkers, 1964; Here, t is time, x is the longitudinal coordinate originating at the landward boundary, h(x, t) is water surface elevation referred to mean sea lev- is the tidal channel width at the water surface ( Figure 2d). A(x, t) is the cross-sectional area of the tidal channel, β cor (x, t) is a correction coefficient arising from the nonuniform distribution of the velocity across the section, r s (x, t) is a term accounting for the water storage over the intertidal platforms adjacent to the tidal channel, and S f (x, t) is the friction slope, that is, the energy dissipation per unit length and unit fluid weight (dimensionless). The two latter quantities read where K T is the conveyance factor, inversely proportional to the flow resistance (Chaudry, 2008;Cunge et al., 1980).
Within the 1D modeling framework adopted here, the continuity Equation 1 accounts for mass conservation within the whole transect area A T (x, y, t). Conversely, the contribution to momentum of the intertidal areas adjacent to the tidal channel is assumed to be negligible. Hence, Equation 2 describes the flow within the active channel area A(x, t). Intertidal areas affect this flow only indirectly, through the correction coefficient β cor and the conveyance factor K T appearing in the friction slope equation. Indeed, in the composite cross sections in the channeled parts of the tidal basin (Figure 2d), the overall flow conveyance (and hence the flow resistance) must account for different values of Gauckler-Strickler coefficient in the main channel (k s ) and over the intertidal areas (k s,l ). Following Engelund ( where A l is the area of lateral intertidal platforms (Figure 2d), R H,l is the corresponding hydraulic radius, while R H is the hydraulic radius of the main channel. Equations 1 and 2 also hold in wide, almost rectangular sections located upstream of the tidal channel head (Figure 2c)   The pressure forces acting at the interface between the main channel and the lateral platforms (black dashed inclined lines in Figure 2d) could affect momentum dissipation when water spills onto the intertidal areas or, vice versa, flows back into the channel. These effects, here neglected as a first approximation, have been investigated by Brown and Davies (2010) with reference to the role of an estuarine channel-sandflat system in producing flood/ebb-dominant net sediment transport and, more recently, by Hepkema et al. (2018), with respect to the sensitivity of tidal bar wavelength to channel width variations.
In the above described 1D modeling framework, the wetting and drying processes typically occurring over the shallower areas of the tidal basin are managed using the physically based procedure proposed by Defina (2000). A subgrid model is introduced to account for ground irregularities and the definitions of conveyance K T , cross-sectional area A, coefficient β cor , and storage coefficient r s are modified to deal with small depth flows (see Table S1 and Figure S1 of SI and details therein). The system formed by Equations 1 and 2 needs to be completed with suitable initial and boundary conditions for the water discharge Q(x, t) and the water surface elevation h(x, t). At the beginning of each simulation, the tidal basin is supposed to be at rest and fully wet.
where h 0 is the still-water level equal to the mean high water level (MHWL). Assuming that the flow is subcritical, a no-flux condition is imposed at the upstream boundary (x = 0) of the tidal basin while downstream, at the channel mouth (x = L), the water level is set equal to the tidal oscillation, namely where a 0 is the tidal amplitude and T is the tidal period.
The initial bathymetry consists of a transversal cosine-like incision of small amplitude, whose minimum elevation is placed along the longitudinal axis of an otherwise flat basin. Some tests on the robustness of the 1D hydrodynamic model are reported in the SI ( Figures S2-S5), where a comparison is carried out with a two-dimensional finite element model (Defina, 2000) which solves the full set of two-dimensional shallow-water equations.

The Cross-Sectional Shear Stress Distribution
The 1D hydrodynamic model provides a cross-sectionally averaged description of the longitudinal bed shear stress. For a given bathymetry, we are however interested in the local value τ(x, y, t) of the stress, which drives bed evolution. This value is here computed through an approach widely adopted in both river (Diplas, 1990;Pizzuto, 1990) and tidal morphodynamics (Fagherazzi & Furbish, 2001). At each transect (i.e., for a given x) and at each tidal-cycle instant (i.e., for a given t), the lateral distribution of longitudinal bed shear stresses is given by the following equation (Lanzoni & D'Alpaos, 2015) where dn is the wetted perimeter between two normals to the bed and dA T is the corresponding area. The dimensional coefficient ( , ) x y  (m 2 ) weights the lateral transport rate of downstream momentum due to the nonuniform along-depth distribution of turbulent stresses. In particular,  turns out to be proportional to the square of the local flow depth, computed along the direction ζ normal to the bed (Figure 2d), and to the local flow conductance, s T k dA dn g (Lanzoni & D'Alpaos, 2015).
Equation 7 shows that the shear-stress distribution across a transect is controlled by two contributions. The first term depends on the cross-sectional geometry, embedded in the term dA T /dn and, indirectly, on the friction slope S f (recall the second of [ with D the local flow depth. Given the water level in a transect, the larger the discharge, the higher the friction slope and, therefore, the bed shear stress. On the other hand, the larger the conductance (and hence the coefficient ), the higher the momentum exchange between the various subareas in which the cross section has been divided. The second term on the right-hand side of Equation 7 is negative within the channel and, therefore, it subtracts momentum from the channel, reducing the shear stresses. In turn, the subtracted momentum concurs to increase the bed shear stress at the edge of the intertidal platforms flanking the tidal channel. Overall, the momentum exchange smooths out the distribution of bed shear stress over the transect (D'Alpaos et al., 2006). In the following, we will use this distribution to determine erosion and deposition processes within any basin transect.

Morphodynamics
The morphodynamic module determines the evolution in time of bed elevations within the tidal basin z b (x, y, t), referred to MSL. Elevations are assumed to change on a temporal scale longer than the tidal period.
In other words, elevations are kept constant throughout a tidal cycle during which the model computes the hydrodynamics, the transverse shear stress distribution, the tidally averaged deposition and erosion rates. Deposition and erosion rates are then used to update the bed introducing a morphological factor, as described below.
The tidal basin is assumed to be composed of fine sediments and, hence, bed load is neglected as compared to suspended load. Moreover, the median sediment diameter is taken to be greater than 20 μm, such that also flocculation can be neglected (Mehta et al., 1989;Winterwerp, 2007). The tidally-averaged bed evolution equation then takes the form (e.g., Marani et al., 2010;) in which ( , , ) x y t  and ( , , ) x y t  are the tidally averaged erosion and deposition rates, respectively, while R is the rate of relative sea level rise (RSLR) accounting for both eustatism and subsidence, assumed constant in space and time.
The presence of vegetation modifies turbulent structures within vegetated patches (Mudd et al., 2010;Vandenbruwaene et al., 2011), changing the intensity of bed shear stresses and, consequently, affecting both erosion and deposition rates. When vegetation encroaches the intertidal platform, the sediment bed is stabilized (i.e., no erosion occurs) and overall sediment deposition is enhanced in relation to the amount of vegetation biomass as will be clarified later. Otherwise, for bare soil, the instantaneous erosion rate is computed through the relation proposed by Carniello et al. (2012) where τ ce is the critical shear stress for erosion, e r0 is an empirical parameter depending on sediment features, while a is a dimensionless parameter. Assuming that τ ce is distributed according to a log-normal probability density function, the previous relation modifies the classic relation proposed by Metha (1984) introducing a smooth transition between τ ≤ τ ce and τ > τ ce . The parameter a controls the degree of smoothing and here is set equal to 2 based on the analyses carried out by Carniello et al. (2012). The tidally-averaged erosion rate  is eventually computed as The overall tidally averaged sediment deposition rate  is, in general, the sum of three contributions where Q s is the deposition rate due to sediment settling, Q t is the deposition rate due to vegetation trapping, and Q o is the organic production rate. Clearly, in the absence of vegetation, only the settling rate contributes to the overall deposition. The settling, trapping, and organic production rates are computed as (e.g., Marani et al., 2010;Temmerman et al., 2003) where ρ s is the sediment density, w s is the settling velocity, Q b0 is the maximum organic production rate, b(z b ) and b max are the local vegetation biomass (see Equation 18) and its maximum value, respectively, α and β are empirical coefficients Mudd et al., 2004), and C is the depth-averaged suspended sediment concentration.
The latter quantity is determined for a typical tidal cycle, by solving the instantaneous sediment balance equation for the water column. Assuming that a given sediment concentration C 0 is supplied externally, this balance equation takes the form Randerson, 1979;Temmerman et al., 2003): where α b β is the trapping increment due to the presence of vegetation, parameterized through its biomass b. The last (source/sink) term on the right-hand side of Equation 16 represents the sediment exchange between the water column and its surroundings (see the detailed derivation in the SI) and reads Randerson, 1979) where dh/dt the rate of change of the free water surface due to the tide. While dh/dt is defined locally at each cross section, the external sediment input C 0 is assumed constant within the basin and equal to the concentration imposed at the channel mouth.
Equation 16 states that the suspended sediment concentration increases in the flood phase and decreases in the ebb phase according to the temporal variation of the local water surface, water depth, and vegetation biomass. Actually, the suspended concentration would change also because of flow advection and diffusion, decreasing progressively in the landward direction and moving away from the main channel. These spatial variations due to advection and diffusion are neglected in the present modeling approach. We will discuss the consequences of this assumption in Sections 3.2, 3.3 and 4.2.

Vegetation Modeling
Halophytic vegetation starts to grow when the elevation of the intertidal platform exceeds a species-dependent threshold. Vegetation biomass generally depends on plant species, soil salinity, hydroperiod, and nutrient availability (Silvestri & Marani, 2004). As a first approximation, biomass can be expressed in terms of platform elevation (Morris et al., 2002). Denoting by b the weight per unit surface of the above-ground vegetation biomass (Mudd et al., 2004), we follow Marani et al. (2013) and write Equations 18 and 19 allow one to describe biomass growth according to the rate λ 2 , its peak at elevation z b0 , and its following decay according to the rate λ 1 , as the platform elevation increases and the hydroperiod decreases. Figure 4 shows an example of three biomass distributions obtained by considering different shapes of the function f(z b ). They can be assumed to represent different vegetation types within the tidal frame. The first vegetation type (type A) is characterized by a rapid biomass growth (λ 2 = 1,000 m −1 ) up to a peak at 0.01 m above MSL (z b0 = 0.01 m above MSL), and then by a slow decay (λ 1 = 15 m −1 ). Biomass vanishes when marsh elevation reaches MHWL. In the second vegetation type (type B), biomass grows at a slightly slower rate before decaying and vanishing at MHWL (λ 1 = 15 m −1 , λ 2 = 80 m −1 , z b0 = 0.01 m above MSL). The last vegetation type (type C) is characterized by a biomass that grows even more slowly and reaches its peak at MHWL, its decay starting for elevations larger than MHWL (λ 1 = 3 m −1 , λ 2 = 18 m −1 , z b0 = 0.50 m above MSL). The first two vegetation types mimic a monospecific vegetation scenario, with halophytes that start growing at their maximum productivity as soon as marsh elevation exceeds MSL, and progressively decrease their biomass as marsh elevation increases toward MHWL (like in Spartina dominated scenarios, Mudd et al., 2004;D'Alpaos et al., 2007). Conversely, the third vegetation type mimics a multispecies vegetation scenario, typical of marshes populated by a variety of vegetation species (e.g., D'Alpaos et al., 2007).
Overall, vegetation affects not only sediment erosion and deposition, as explained above, but it also increases flow resistance, favoring flow concentration within the channels (Temmerman et al., 2007). Moreover, when vegetation starts encroaching the tidal platform, a logarithmic velocity profile does not hold anymore (Nepf, 1999;Nepf & Vivoni, 2000) and the estimate of the transverse exchange of momentum in Equation 7 becomes problematic. We then directly embed this effect in the Gauckler-Strickler coefficient (recall that  s k  ) which is decreased to account for the drag force due to plant stems. The overall Gauckler-Strickler coefficient k s on the vegetated platform is hence computed as (Belliard et al., 2015): where k s,g accounts for grain roughness while k s,v quantifies the decrease in conductance (i.e., increase in friction) due to vegetation. No bedform roughness is considered, owing to the fine, almost cohesive character of the sediment. The larger vegetation biomass, the larger the friction, that is, the smaller k s . The coefficient k s,v is thus mediated by the ratio b/b max between the local and the maximum plant biomass.

Three-Dimensional Domain Schematization
The tidal basin geometry is described through M transects with shape given by the bed elevation z b (x, y). In order to study tide propagation over the basin using a 1D model, each cross section is approximated by using the compact/compound framework described in Section 2.1.
Typically, the tidal channel region occupies the portion of the basin where the streamlines are roughly straight. Various criteria may be used to separate the tidal channel area from the lateral flow regions (Fagherazzi et al., 1999;Geng et al., 2018). As previously pointed out, we assume that intertidal areas adjacent to the tidal channel concur only to water storage in the continuity equation. We thus use a simple geometrically based criterion to identify the tidal channel region. Specifically, a given bed point of a transect is assumed to belong to the tidal channel when the angle between the normal to the bed and the vertical exceeds a given threshold (say 1°), that is when the local transversal slope exceeds 1%-2%.
After identifying the unchanneled and channeled regions within the basin, for each transect we compute the channel width ( , ) x t  and the mean elevations of the channel z b,ch (x, t) and of the lateral areas z b,l (x, t). Next, we assign the local values of the Gauckler-Strickler friction coefficients k s (x, y, t), according to Equation 20.

Design of Numerical Experiments
Each simulation proceeds as follows. We first compute the 1D flow field that establishes in the basin under the action of the seaward imposed tidal forcing. At the other three boundaries, no flux conditions are imposed and, during the whole tidal cycle, the bed is kept fixed. The 1D hydrodynamic Equations 1 and 2 are solved numerically by the Preissman method (Cunge et al., 1980) using a timestep of 200 s. For a given bed topography, two tidal cycles are simulated to ensure a full adaptation of the flow field to the bed configuration. We then use the values S f (x, t) of the friction slope provided by the 1D model to determine the cross-sectional shear stress distribution at each instant of the tidal cycle. To this aim, Equation 7 is solved numerically by centered finite differences (Lanzoni & D'Alpaos, 2015;Pizzuto, 1990). Next, we use the water column Equation 16 to compute the tidally averaged values of the local erosion and deposition rates. When the platform is vegetated, vegetation increases friction in the hydrodynamic module (Equation 20), changing the shear stress distribution in Equation 7 through , and affects depositional processes in the morphodynamic module (Equations 14 and 15).
Finally, relying on an offline technique (Roelvink, 2006), the tidally averaged values of erosion and deposition rates are used to update the local bed elevations using a morphological acceleration factor of 30. Erosion and deposition rates after 1 day of simulation are thus used to approximate their values for the next 30 days of simulation. The above-described procedure is repeated until the system reaches an equilibrium configuration. The overall duration of the simulation is set by a test/trial procedure that allows to check whether the basin asymptotically reaches constant bed elevations everywhere.
The input parameters used by default in the model are reported in Table 1. Three different scenarios have been investigated. First, we considered purely erosional conditions. Second, we accounted for erosion, sea level rise, and settling deposition. The latter process strives to balance out the former two, eventually SGARABOTTO ET AL.

Purely Erosive Equilibrium
Model robustness was tested for the case of purely erosional conditions. A sensitivity analysis was carried out varying tidal basin dimensions (L, B T ), initial tidal-flat elevation (z in,tf ), tidal amplitude (a 0 ), critical shear stress for erosion (τ ce ), Gauckler-Strickler coefficient (k s ), and the basin aspect ratio, as in Lanzoni and D'Alpaos (2015) ( Table 2). The simulations were carried out for 1,000 years.
The effects of changes in a given parameter, keeping all the others unchanged, are summarized as follows. An increase in basin width from 100 to 500 m makes the channel deeper keeping the tidal channel width approximately constant. Moreover, the channel deepens and widens by changing the tidal amplitude from 0.1 to 0.7 m. Differently, when the initial elevation of the basin increases from −1.0 to −0.1 m above MSL, the channel becomes narrower and shallower. Similarly, increasing the critical shear stress for erosion from 0.2 to 0.7 Pa makes the channel shallower and narrower. Figure 5a shows an overall view of the spatial distribution of bed elevations, the distribution of longitudinal maximum flow velocity and longitudinal maximum bed shear stress computed for the final equilibrium configuration of a typical test (I3 in Table 2), with initial conditions depicted in Figures 5b-5d. This test is taken as reference for the comparisons carried out in Figures 5e-5n, reporting the along-channel width distributions and the longitudinal channel bed profiles computed in the other tests of Table 2. Overall, width distributions show an almost linear trend, except close to the channel head where the bed is subjected to wetting and drying during the tidal cycle. There, width decreases faster than linearly. The longitudinal profile shows an approximately linear trend except toward the channel head, where it joins the landward intertidal platform with a sharp transition.
As the tidal basin widens, the main channel adapts to the larger tidal prism by becoming deeper while maintaining approximately the same width distribution (Figures 5e and 5f). Moreover, the larger the tidal amplitude, the deeper and wider the channel (Figures 5g and 5h). For tidal amplitudes larger than that of the reference test (0.5 m) the channel becomes larger and deeper (Figures 5g and 5h). Moreover, the channel lengthens through a landward migration of the step-like transition at the channel head. The opposite behavior is observed for tidal amplitudes smaller than 0.5 m.
The higher the initial elevation of the tidal flats and, hence, the shallower the initial configuration of the tidal basin, the shallower the channel (Figures 5i and 5j). Moreover, increasing the initial elevation of the tidal flats, the channel shortens because of the seaward migration of the channel head. The trends observed for the channel bed profile for deeper initial basins are similar to those observed for wider tidal basins. As this critical threshold increases, erosion decreases making the channel narrower, shallower and shorter owing to the seaward retreat of the channel head (Figures 5k and 5l). Note that, for a critical threshold higher than SGARABOTTO ET AL.   Table 2). Arrows highlight morphodynamic tendencies. The ranges of investigated parameters are those reported in Table 2. The basin width T  ranges between 100 and 500 m, increasing by 100 m in each test. Tidal amplitude a 0 changes between 0.3 and 0.7 m, with an increase of 0.1 m in every test. Critical shear stress τ ce ranges between 0.2 and 0.6 Pa, increasing by 0.1 Pa in each test. Initial bed elevation z tf,in changes between −1.0 m above MSL and −0.1 m above MSL with an average increase of 0.2 m in every test. Gauckler-Strickler coefficient ranges between 20 and 40 m 1/3 /s, with an increase of 5 m 1/3 /s in every test. that of the reference case, downstream of the channel head, the bed profile keeps an almost linear trend with minor changes because of the reduced erosion rates. Qualitatively similar trends are observed as the Gauckler-Strickler resistance coefficient increases (Figures 5m and 5n).
Finally, the larger the Gauckler-Strickler coefficient, the shallower, shorter, and narrower the channel (Figures 5m and 5n). Indeed, as the Gauckler-Strickler coefficient increases, the friction decreases and, hence, erosion reduces.

Morphodynamic Equilibrium in the Presence of the Competing Effects of Deposition, Erosion and Sea Level Rise
The present model, besides erosional effects, accounts also for the effects of RSLR and deposition of suspended sediment, as described in Section 2.3. We are interested in assessing under which conditions a dynamic equilibrium can be attained both on the tidal channel and the adjacent intertidal platforms, depending on the rate of RSLR and the external availability of suspended sediment, embodied by the prescribed concentration C 0 . In order to better clarify the role played by the various mechanisms, here we focus only on settling deposition and postpone to the next section the analysis of vegetation effects. In addition, for comparison purposes, we refer to the parameter set of run I2 in Table 2. However, here we consider deposition effects and assume a rate of RSLR equal to 3.5 mm/yr, which is a reasonable value for the Venice lagoon .
The results obtained by supplying externally a concentration C 0 = 30.0 mg/l are first considered. The longitudinal profiles in Figure 6a indicate that, in the early stages of evolution (after about 50 years), a steplike transition forms at the conjunction between the landward platform and the channel head. There, the change in active cross-sectional area leads to an almost flat segment break in the longitudinal distribution of the tidal prism, which otherwise would increase linearly seawards (Figure 6b). As the landward platform SGARABOTTO ET AL.

10.1029/2020WR028577
13 of 24  (Figures 6c and 6d), the tidal prism continuously decreases in time. After about 500 years, the depositional and erosional processes offset each other, the tidal prism distribution sets at approximately steady values. The active channel bed profile tends to assume a concave shape while a convex shape characterizes the bed profile upstream of the active channel head (Figures 6a and 6d). The active channel bed deepens up to about 100 years and then aggrades until it reaches the final equilibrium condition. Indeed, as the tidal prism decreases in time, the velocity within the active channel also decreases. Consequently, deposition starts to prevail over erosion and bed aggradation occurs. Furthermore, since the beginning of the evolution process, deposition prevails in the landward portions of the basin and outside of the active channel leading to higher bed elevations. This is due to the quite slow water currents that characterize those areas, where erosion is almost negligible. The bed elevation of tidal flats adjacent to the active channel increases faster as moving seawards. Moreover, the upstream tidal flat merging with the channel head exhibits negligibly small variations of bed elevation in the transverse direction (Figures 6c and 6d).
The route toward an equilibrium condition is evident looking at the temporal evolution of mean elevations of the channelized region (CH) and the tidal flat region (TF) (Figure 7). The consequences of varying the sediment input are also shown in Figure 7. For a sediment input of 30.0 mg/l, erosional effects are found to prevail at the beginning of the simulation. The initially channelized region progressively deepens, stretching throughout the basin, and assumes a funnel shape (Figure 7c). Meanwhile, the adjacent intertidal platform, differently from the purely erosive setting, gradually increases its elevation due to sediment settling. This accretion is also related to the small values attained by bed shear stresses on intertidal platforms, implying a negligible erosion rate. Within the channel region, where shear stresses and, hence, erosional effects are stronger, deposition and erosion processes slowly offset each other. Bed variations throughout the basin become progressively smaller and the system finally tends toward an equilibrium configuration. Intertidal platforms flanking the channel reach a mean elevation of 0.16 m above MSL, while the channel region converges to an elevation of −2.49 m above MSL. The lateral tidal flats reach their equilibrium elevation (0.16 m above MSL) around 400 years before the channel bed (−2.49 m above MSL). Due to the concurrent presence of both depositional and erosional effects, the final equilibrium is dynamic, that is, erosion balances deposition. The basin infilling produces a narrowing of the channel (Figure 7b). After the rapid changes at the beginning of the simulation, the mean width decreases until it reaches its minimum value (around 19 m) at equilibrium.  increases (C 0 = 45.0 mg/l) the equilibrium elevation of the intertidal platform increases to 0.27 m above MSL, and the equilibrium depth of the channel area reduces to −2.1 m above MSL, due to the decrease in the tidal prism. On the contrary, a decrease in sediment availability (C 0 = 15.0 mg/l) leads to a decrease in platform elevation below MSL (−0.2 m above MSL) and a deepening of the channel (−3.1 m above MSL) due to the increase in the tidal prism. Interestingly, the time required to reach the equilibrium configurations decreases as sediment availability increases (Figures 7a and 7b), that is, systems with larger sediment supply equilibrate faster to changes in the environmental forcing than systems with lower sediment supply (e.g., D'Alpaos et al., 2011). In particular, the delay between the time required by the active channel and the tidal flats to reach the equilibrium configuration decreases from around 100 years for C 0 = 15 mg/l to ∼350 years for C 0 = 45 mg/l.

Ecomorphodynamic Equilibrium: The Role of Vegetation Growth
We now consider the effect of vegetation, with reference to the biomass distribution of type A shown in Figure 4. The other relevant parameters are those of Test Case I2, R = 3.5 mm/yr and C 0 = 30.0 mg/l. The longitudinal distributions of channel thalweg, local elevation of marsh edges, and tidal prism obtained in the presence of type A vegetation are shown in Figure 8.
The results indicate an evolution trend qualitatively similar to that obtained in the absence of vegetation ( Figure 6). A step-like transition is again observed to form at the channel head, followed by a weakly convex up thalweg profile. No differences occur as compared to the unvegetated scenario until the bed elevation is below the threshold for vegetation growth. This threshold is first exceeded by the landward intertidal platform, which accretes faster owing to the larger local deposition rate (see  platform and the tidal channel (Figure 9). In particular, it leads to a faster reduction in the tidal prism. As a consequence, the discharge flowing through the tidal channel decreases, contributing to channel infilling. However, vegetation growth also increases flow resistance on the lateral platforms and promotes flux concentration and erosion within the channel. Vegetation encroachment enhances marsh accretion through sediment trapping and organic soil production. In the case of C 0 = 30.0 mg/l, this boosted accretion allows the marsh to reach a maximum elevation of ∼0.19 m above MSL after about 300 years (Figure 9a). A faster platform accretion, in turn, feeds back into channel bed evolution. First, the main channel becomes deeper. However, when sediment deposition becomes dominant, the channel bed accretes reaching a mean elevation of ∼−1.80 m at 250 years. Then, vegetation favors flux concentration within the channel due to increased friction on the tidal platform adjacent to the channel. Erosion increases again, offsetting depositional effects and leading asymptotically to a mean channel bed elevation of ∼−2.64 m. The channel reaches equilibrium about 600 years later than the marsh platform. Furthermore, on average, the final channelized region is not only deeper but also narrower (16 m) as compared to the correspondent unvegetated case (19 m).
In the case of a larger sediment supply (C 0 = 45.0 mg/l), the intertidal platform further increases its elevation and the channelized region becomes on average shallower than in the unvegetated case. However, the channel is as wide as in the unvegetated case because the flux is already mainly confined within the channel when vegetation starts to grow. On average, channel equilibrium and marsh platform equilibrium are delayed by ∼400 years. Conversely, in the case of a lower sediment supply (C 0 = 15.0 mg/l) the platform does not emerge, and vegetation cannot encroach the platform. Hence, the channel bed tends to the same equilibrium configuration as in the correspondent unvegetated case. In this case, the channel equilibrium is delayed by ∼50-100 years with respect to marsh platform equilibrium.
A direct comparison between vegetated (type A) and unvegetated scenarios is pursued in Figure 10, which shows the evolution in time of the mean elevations within the channel and on the intertidal areas. For low sediment supply (C 0 = 15.0 mg/l), the mean elevations of the channel region and of the marsh region are the same for both scenarios (Figure 10a). Indeed, the mean elevation of the intertidal platforms keeps below the threshold at which vegetation starts to grow. When sediment supply is increased to 30.0 mg/l, the intertidal platforms emerge above MSL, promoting vegetation growth. This enhances platform accretion, which is slightly higher than in the unvegetated case. In turn, the increased friction on the vegetated platform concentrates the water flux into the channel and, hence, enhances channel deepening (Figure 10b). A further increase in the sediment supply (e.g., C 0 = 45.0 mg/l) boosts intertidal platform elevation, whereas vegetation growth promotes channel deepening compared to the unvegetated case. Nevertheless, when the SGARABOTTO ET AL. sediment supply is relatively high, vegetation mildly affects the equilibrium platform elevation, which is mainly controlled by settling deposition. The deepening of the channel is thus due to flux concentration induced by increased friction on vegetated areas (Figure 10c).
The role of the different contributions to marsh accretion (i.e., settling, trapping, and organic production) expressed as mm/yr, is presented in Figure 11. For low sediment supply (C 0 = 15.0 mg/l), settling deposition is the only contribution to total accretion because intertidal-platform elevation does not exceed the threshold elevation for vegetation growth (Figures 11a). When sediment supply increases (C 0 = 30.0 mg/l), the intertidal platform reaches an elevation that allows vegetation growth. Organic production and trapping deposition seemingly equal their intensity, while the settling contribution is four times larger than the other two contributions (Figures 11b). When the averaged biomass growth reaches its maximum within the basin, the total accretion rate reaches a maximum too. The effects of vegetation become milder for higher sediment supply (C 0 = 45.0 mg/l, see Figures 11c). Indeed, in the latter case, sediment settling largely overwhelms trapping deposition and organic production and largely contributes to the final elevation of the marsh platform and of the related channel geometry. Note that organic production and trapping deposition SGARABOTTO ET AL.  peak at mean elevations lower than MSL because of the longitudinal gradients in marsh-platform elevations (Figures 8c and 8d). As previously discussed, vegetation first colonizes the landward portions of the basin, when the mean elevation of the intertidal platform is still below the threshold for vegetation growth. In the following stages, seaward intertidal areas progressively reach the threshold elevation for vegetation growth. This leads to a progressive increase in the mean elevation of the vegetated areas, as shown in Figure 10.
The evolutionary scenario described by Figures 8 and 9, as well the related vegetation feedback mechanics, can obviously be modified if another type of vegetation (e.g., type B and C) is considered. We will discuss in the next section how different types of vegetation may affect the width-to-depth ratio of the channel cross sections.  (Marani et al., 2002). Panel (a) shows the relationship between channel cross-sectional area and tidal prism. The continuous line represents the slope (6/7) of the O'Brien-Jarrett-Marchi law (D'Alpaos et al., 2009). Panel (b) represents the relationship between the width-to-depth ratio and channel width. Field data from Marani et al. (2002) are denoted by crosses if referred to salt-marsh channels, and by plus symbols if referred to tidal-flat channels. The corresponding fitting lines are denoted in gray (MA) and black (TF), respectively. Panel (c) shows box plots of the width-to-depth ratio obtained with the present model and observed in the field. Black circles in panels (a) and (b) represent the results obtained considering the purely erosional setup; red circles represent the results obtained considering the depositional model; blue circles represent the results obtained considering the vegetation model. settling deposition, vegetation effects (trapping, organic soil production, increase in bottom roughness over the vegetated surface) and RSLR (vegetated scenario). The comparison is carried out in terms of: the relationship between channel cross-sectional area, A, and the tidal prism, P (Figure 12a); the relationship between the width-to-depth ratio, β, and channel width, B (Figure 12b); and box plots of β values (Figure 12c). Field data referring to tidal channels in the Venice Lagoon (Marani et al., 2002) are also shown in Figures 12b and 12c. Model results in the case of purely erosive conditions are quite similar to those obtained through LD15's model (Lanzoni & D'Alpaos, 2015). Indeed, the simulated tidal basin is short enough to ensure that the assumption of a quasistatic propagation of the tidal wave holds and, hence, momentum exchange between the various channel sections, accounted for in the present model, are likely to play a minor role. In particular, the A-P relationship follows a classical power law of the type A ∝ P α with α = 6/7 (D'Alpaos et al., 2009) throughout the entire basin (Figure 12a, black circles), suggesting the reliability of the computed channel sizes in response to the landscape-forming tidal fluxes Friedrichs, 1995;van der Wegen et al., 2010).

Overall Comparison of the Results
When the effects of sediment supply, sea level rise, and vegetation growth are taken into account, a power-law relationship still holds between A and P. However, deviations clearly emerge for the landward cross sections, far from the inlet and highly subjected to wetting and drying processes. Model results support previous speculations put forward in the context of modeling frameworks which did not account either for wetting and drying processes or vegetation effects on channel cross-sectional shape and size Lanzoni & D'Alpaos, 2015). The relationship between β and B (Figure 12b), shows that, under purely erosive conditions (black circles), the synthetic channel cross sections are wider than those in the other considered cases and are characterized by larger width-to-depth ratios (10 < β < 45). The latter values of β are in good agreement with those observed by Marani et al. (2002) for the tidal-flat channels in the Venice Lagoon (8 < β < 50). Conversely, when the effects of sediment deposition (red circles) and vegetation growth (blue circles) are accounted for, channel cross sections become narrower and deeper (6 < β < 10 and 5 < β < 10 for the depositional and vegetated scenarios, respectively). In this case, β-values compare well with those observed by Marani et al. (2002) for the salt-marsh channels of the Venice Lagoon (5 < β < 8).
In addition, the box plots of Figure 12c show that, for purely erosive conditions (P in Figure 12c), the synthetic cross sections exhibit a median value of the width-to-depth ratio (14.3) which closely approximates that characterizing the tidal-flat channels (CTF in Figure 12c) in the Venice Lagoon (15.1). However, the outliers of the box plots indicate that the variety of channel morphologies observed in the tidal-flat channels of the Venice Lagoon is larger than that computed by the model for purely erosive conditions. On the other hand, in the depositional and vegetated scenarios, the synthetic width-to-depth values better describe the channels cutting through the salt marshes of the Venice lagoon. Indeed, not only the observed median value of β (6.0) is quite close to that characterizing the synthetic morphologies (6.5 and 6.9 for the vegetation and depositional scenarios, respectively), but also the observed and computed variances are very similar (1.9 for the observed β-variances as compared to 0.9 and 0.8 for the depositional and vegetation models, respectively).
The important role exerted by vegetation on the width-to-depth ratio clearly emerges from the box plots of Figure 13, which shows β-values for the three different vegetation types considered in our analyses (see Figure 4). For low sediment supply (C 0 = 15 mg/l), the depositional and vegetated scenarios lead to equivalent results (median β values all equal to 7.3 Figure 13a). Indeed, in this case, the low intertidal platform elevations prevent vegetation growth and the plano-altimetric equilibrium configurations of the different cases are the same. When sediment supply increases (C 0 = 30 mg/l, Figure 13b), the monospecific vegetation types (types A and B) encroach the platform as soon as local elevation exceeds MSL and rapidly increase their biomass. Vegetation influences the flow field, promoting flux concentration within the channels, which become deeper than in the (unvegetated) depositional scenario. A median width-to-depth ratio of 6.4 characterizes channel cross sections for both vegetation types. This value of β is quite similar to that of the median width-to-depth ratio observed for salt-marsh channels in Venice (6.0) and smaller than the value (6.9) which characterizes the purely depositional scenario. Conversely, in the multispecies vegetation case (type C) biomass production is large only when the platform reaches elevations close to MHWL. In this case, the reduction in the tidal prism due to high marsh elevations, leading to channel infilling, prevails over the within-channel flux-concentration effect due to vegetation growth, and the channels display larger width-to-depth ratios. Overall, in the case of type-C vegetation, channel cross sections display a median β-value of 6.8, quite close to that obtained for the purely depositional model (6.9). When sediment supply further increases (C 0 = 45.0 mg/l, Figure 13c) in all considered vegetated and unvegetated cases, marsh equilibrium elevation is so high in the tidal frame that the effects related to the reduction in the tidal prism prevail over vegetation effects. The width-to-depth ratio is equal to 6.2 in the unvegetated and type A vegetation case; to 6.1 in the type B vegetation case; to 6.5 in the type-C vegetation case. In conclusion, vegetation produces two counteracting effects, which are both captured by the proposed framework, differently from previous studies (e.g., D'Alpaos et al., 2006;Temmerman et al., 2007). On the one hand, vegetation tends to concentrate tidal flows within the channel, promoting channel incision (Temmerman et al., 2007) and the formation of deeper channels with smaller width-to-depth ratios. On the other hand, vegetation enhances marsh accretion and the reduction in the tidal prism, fostering shallower cross sections with larger widthto-depth ratios (D'Alpaos et al., 2006). These two counteracting effects may rule out one another, depending on the elevation of the marsh platform flanking the channels and on biomass distribution in the tidal frame. For large biomass productions at relatively low marsh elevations (see Figures 4a and 4b), vegetation effects on within-channel flux concentration prevail in the case of marshes slightly higher than MSL. Conversely, the reduction in the tidal prism dominates in the case of marshes at high elevations in the tidal frame. Furthermore, for large biomass productions at relatively high marsh elevations (see Figure 4c) vegetation-induced flux concentration slightly influences channel cross-sectional size and shape.

Model Limitations
The proposed modeling approach reasonably predicts the morphology of salt-marsh channels and relatively short (i.e., a few kilometers long) tidal-flat channels. The model is deemed to capture the mutual biogeomorphic interactions between tidal channels and the adjacent intertidal platforms. Clearly, tidal channel morphologies can be only qualitatively predicted when model limits of applicability are not strictly met. We briefly recall here these limitations.
Analogously to exploratory models (Murray et al., 2008;Townend, 2010), the proposed framework allows an efficient and thorough exploration of the parameter space, in terms of tidal forcing, external sediment supply, rate of sea level rise, and different vegetation types. Despite the simplifications introduced with respect to morphodynamic models solving the full set of two-dimensional equations Dam et al., 2016;Hibma et al., 2003;, the present model proves quite robust when compared to field observations. Finally, the model allows one to unravel the effects of vegetation growth on the equilibrium configuration of tidal channels cutting through salt-marsh platforms.
SGARABOTTO ET AL. In the box plots "D" refers to the depositional model; "V" refers to the vegetation models in the type A, B, C scenarios; "CTF" refers to tidal-flat channels and "CMA" to salt-marsh channels in the Venice lagoon (Marani et al., 2002).
As to model limitations, we first notice that streamlines over the actual intertidal areas are not parallel to the longitudinal axis for some instants of the tidal cycle. Therefore, the hydrodynamic module cannot simulate streamline deflection toward the channel during the initial stages of the flood phase and the final stages of the ebb phase (Lawrence et al., 2004). Thus, the model cannot simulate either the development of branching channel networks or channel competition for draining the tidal landscape (van Manen et al., 2013;Xu et al., 2017).
Advection of suspended sediment and spatial settling lag effects are ignored in the concentration Equation 16. Furthermore, the external sediment supply during the flood phase, C 0 , is assumed constant within the tidal basin. Consequently, the tidally averaged sediment concentration might be overestimated in the landward channel cross sections, where deposition tends to out-balance erosion (Figures 6a and 8a). Similarly, deposition rates are likely to be overestimated toward the lateral boundaries of the tidal basin (Figures 6c, 6d, 8c and 8d). As a result, spatial variations of intertidal platform topography are somewhat smoothed out. Clearly, the shorter the channel, the milder the effects exerted by advection on the equilibrium profile. These effects are likely to become more important for long channels. Note that the importance of internally generated advective sediment transport in tidal embayments (i.e., in the absence of externally prescribed overtides) has been demonstrated to become large especially for tidal embayments with landward diverging widths (Meerman et al., 2019).
Clearly, further processes as wind waves and overtides, not accounted for in the present model, may affect the modeled equilibrium morphologies and merit attention. Wind waves in general enhance erosional processes on the tidal platform  and at the marsh edge (Bendoni et al., 2016;Leonardi & Fagherazzi, 2014;Leonardi et al., 2016;Marani et al., 2010). In turn, the larger tidal prism due to an enhanced platform erosion together with marsh edge erosion, promote channel widening and deepening. Vegetation cover makes salt marshes more resistant to the erosion induced by wind waves (Carniello et al., 2013;Möller et al., 1999Möller et al., , 2014 and, consequently, tends to counteract the erosion of channel edges. The tidal asymmetry induced by externally imposed overtides may determine a net sediment transport which influences the morphological equilibrium. For tidal channels dissecting coastal lagoons, overtides have been demonstrated to provide some corrections to the equilibrium bed profile  and to determine the possible existence of multiple equilibria in fairly long tidal channels (Schuttelaars & de Swart, 2000). For the considered tidal-basin geometry, the effect of tidal asymmetries caused by external overtides is found to produce a wider and shallower tidal channel in the purely erosional scenario ( Figure S7a of SI). This result is in line with the finding of Toffolon and Lanzoni (2010), obtained by using a linerized perturbative framework and indicating a tendency of the equilibrium bed level to increase in elevation when the phase difference between the semi-diurnal tidal constituent M2 and the M4 overtide falls in the rage (−π/2, π/2). Conversely, the M4 overtide is found to produce minor differences in the equilibrium channel configuration when considering the erosion-settling deposition and the vegetated scenarios (Figures S7a and S7b of SI).
Within its limits of validity, the simplified model can be used to predict the long-term coupled biomorphodynamic evolution of salt-marsh platforms and tidal channels cutting through them associated with changes in external forcing. It is worthwhile noting that, in the long-term, process-based models solving the full set of governing equations (Belliard et al., 2015;Hibma et al., 2003;Xu et al., 2017) not always provide reliable information on the evolution of the systems they describe. Indeed, process-based models have been demonstrated to perform well when applied under constant forcing conditions and in a constrained system that limits the possible morphological changes (e.g., Dam et al., 2016). Conversely, results may drift away from observations over long time scales when considering less constrained environments. As an example, Coco et al. (2013) showed that simulations carried out under the same conditions, but with different numerical models, lead to different final configurations.

Conclusions
The present contribution aimed to study the plano-altimetric biomorphodynamic evolution of a tidal channel and the adjacent intertidal areas. This has been investigated by considering three different scenarios, which progressively include the relevant processes. In the purely erosional scenario, only sediment erosion is assumed to shape the channel and the intertidal platform. In the depositional scenario, the effects of sea level rise and settling deposition, driven by externally supplied sediment, are also accounted for. Therefore, the model describes the intertwined vertical accretion of the intertidal platforms and the related dynamic evolution of the channel which dissects them. Finally, in the vegetated scenario, increased flow resistance induced by vegetation, sediment trapping, and organic soil production sum to the previously described processes, dictating the biomorphodynamic evolution of the system.
The main conclusions of our analyses can be summarized as follows. Under purely erosive conditions, the model produces a larger variety of synthetic channel morphologies which, overall, well agree with the morphological features of tidal-flat channels in the Venice Lagoon. Model results reproduce classical relationships of geomorphic relevance, such as the tidal-prism channel-area relationship and the relationship between the width-to-depth ratio, β, and channel width. Modeled β values (10 < β < 45) are in good agreement with those observed by Marani et al. (2002) for the tidal-flat channels in the Venice Lagoon (8 < β < 50). Moreover, the inclusion of wetting and drying effects, as well as the connectivity of the various transects ensured by the newly developed 1D hydrodynamic model, allow one to reproduce the step-like transition that forms at the conjunction between the landward intertidal platform and the channel head (Figure 1).
When the effects of sediment deposition, sea level rise, and vegetation growth are accounted for, model results better agree with the morphological features of observed salt-marsh channels. In particular, the synthetic β values obtained when considering vegetation effects (5 < β < 10) compare quite well with observed ones (5 < β < 8).
Vegetation growth on the marsh surface produces two counteracting effects. First, vegetation increases flow resistance over the marsh platform, thus promoting flow concentration within the channel. This favors channel incision and the development of narrower and deeper cross sections (with lower width-to-depth ratios). Second, the enhanced inorganic deposition and the organic soil production induced by vegetation lead to a faster vertical accretion and to a faster reduction in the tidal prism. This favors channel infilling and the development of shallower cross sections (with larger width-to-depth ratios). Whether the first effect prevails over the second or vice versa, depends on the elevation of the marsh platform in the tidal frame. This elevation is in turn dictated by sediment availability, rate of sea level rise, and vegetation biomass distribution within the tidal frame. For large biomass productions at relatively low marsh elevations (slightly higher than MSL), the flux concentration induced by vegetation prevails. Conversely, the reduction in the tidal prism dominates for marshes in which biomass production is large at high elevations in the tidal frame (around MHWL). In this type of marshes, vegetation-induced flux concentration within the channel slightly affects channel cross-sectional size and shape.