Our starting point is a two-dimensional shallow water model that solves the hydromorphodynamic problem by integrating numerically the depth-averaged shallow water equations coupled with the Exner equation, which describes the time evolution of riverbed elevation41. When paired with a description of the vegetation dynamics and its effect on flow and sediment transport, such a model was demonstrated to reproduce key ecomorphodynamic processes34. In this study, we model feedbacks between hydromorphodynamic processes and vegetation including a description of both above- and below-ground plant traits and their dynamics. In particular, we consider that vegetation impacts flow and sediment transport by modifying the flow resistance, the bed shear stresses, and the threshold for the onset of bed load transport. In turn, morphodynamic processes occurring during flood events are responsible for two main vegetation mortality mechanisms, which are, plant uprooting and sediment burial, while water level fluctuations during low flow periods between floods control the vegetation growth.
Hydromorphodynamic processes
River hydromorphodynamics is simulated using the numerical model BASEMENT (freeware software)41 and the computational domain is discretised using a triangular unstructured grid. First, the model solves the hydrodynamic problem by using the Manning-Strickler approach for evaluating the hydraulic roughness, in which the bottom shear stress, τ, reads as
$$begin{aligned} mathbf{{tau }} = frac{rho g mathbf {u}|mathbf {u}|}{K_s^2h^{1/3}}quad , end{aligned}$$
(1)
where (rho) is the water density, g is the gravitational acceleration, (mathbf {u}) is the flow velocity vector, h the water depth and (K_s) the Strickler coefficient (the inverse of the Manning coefficient n). Second, the sediment continuity equation (Exner equation) is solved to obtain the evolution of the bottom elevation (z_b) of a riverbed composed of a uniform sediment. The bed load transport intensity is evaluated using the standard Meyer-Peter and Müller formula42, as a function of the excess of the Shields shear stress (theta) above a threshold value ((theta _{cr}=0.047)), where
$$begin{aligned} theta = frac{|mathbf{{tau }}|}{(rho _s -rho ) g d_s} end{aligned}$$
(2)
and (rho _s) and (d_s) are the sediment density and diameter, respectively.
Vegetation description
Vegetation is described by a total dimensionless biomass density, B, which is partitioned into two components, above-ground, (B_c) (subscript c stands for canopy) and below-ground, (B_r) (subscript r stands for roots) (Fig. 2a). The above-ground biomass is considered to be evenly distributed along the plant height, H. Conversely, the below-ground biomass is characterised by a dimensionless vertical density distribution (b_r(zeta )) that extends downward in the (zeta)-direction, from the riverbed surface ((zeta =0)) to the rooting depth (zeta _r(t)), i.e. the maximum depth reached by roots at a generic time, t (Fig. 2a). (b_r(zeta )) is calculated via the stochastic model developed by Tron et al.43. The model assumes that fluctuations of the water table level follow the water level in the channel, which produce an alternating sequence of root growth and decay periods at each riverbed depth (zeta). This behaviour is then described as a stochastic process and solved to obtain a steady-state solution for (b_r(zeta )) (see Appendix for the formulation43), which depends on the mean frequency, magnitude, and decay rate of water table fluctuations and represents a key component in our model, controlling plant growth rate and vegetation resistance to uprooting.
Vegetation growth dynamics
The fluctuation of the water table level is one of the main driving factors controlling growth performances of riparian plants22. The ability of species to rapidly grow roots tracking the water table was found to be key for determining the growth rate of plants and their potential establishment success on bars44. Highly variable water table levels during growth periods tend to produce water stress that reduces plant growth because of reduced root respiration, while long free-inundation periods may be similarly harmful for plants that fail to grow roots deep enough to secure a connection with the groundwater. To model this link, we use the function (b_r) as a proxy for the ability of the plant to tolerate inundations, which depends on the riverbed depth reached by the roots a certain time, and relate that to the plant growth rate. We consider that the total biomass density B (above- and below-ground components) grows in time (t) following a logistic function:
$$begin{aligned} frac{dB}{dt} = sigma _B B biggl (1-frac{B}{B_{max}}biggr ) quad , end{aligned}$$
(3)
where (B_{max}) is the maximum biomass value (set to 1 in our model) and (sigma _B) is the growth rate assumed to vary depending on the dimensionless root density distribution, (b_r), as
$$begin{aligned} sigma _B = phi _{B} int _{0}^{1}b_{r}(z)dz end{aligned}$$
(4)
where (phi _{B} [-]) is a scaling factor and (int _{0}^{1}b_{r}(z) dz) represents the steady-state dimensionless root biomass. With Eq. (4), we assume that vegetation grows faster (higher (sigma _B)) when plant roots are more developed (at the steady-state). This assumption is largely used in modelling dryland vegetation where plant species act as phreatophytes45, similarly to riparian plants. Biomass density decay due to waterlogging is included considering that B decreases exponentially with a constant rate of 0.1 when the riverbed level falls below the mean water table level.
We assume that the rooting depth changes over time following an exponential function as
$$begin{aligned} frac{dzeta _r}{dt} = sigma _{r}( zeta _{r,max}-zeta _r) quad , end{aligned}$$
(5)
where (sigma _r) is the root deepening rate, which is constant, and (zeta _{r,max}) represents the distance between the riverbed surface and the minimum water table level (Fig. 2a). Below such level, riverbed matrix is saturated with pore water and roots cannot grow due to the resulting anoxic conditions43. Equation (5) implies that roots grow faster as farther they are from the groundwater and linearly with (zeta _{r,max}). This behaviour is representative of phreatophytes plant species that uses groundwater as main source of water and tend to elongate roots to keep pace with the receding rate of the water table level21,46.
The proportion of the total biomass growth allocated to each plant component is derived using a mass balance model47 by introducing two constant partitioning coefficients, (lambda _i) with (i in (c,r)) (i.e. canopy and roots). These define the fraction of the total biomass growth allocated above- and below-ground and satisfy (lambda _c + lambda _r = 1) for all times. As a consequence, the growth rates of the two plant components can be written as
$$begin{aligned} frac{dB_i}{dt} = lambda _i frac{dB}{dt} ;. end{aligned}$$
(6)
Here we consider (lambda _i) constant, assuming no plasticity in biomass partitioning. The canopy height depends on the above-ground biomass through the function48
$$begin{aligned} H(t)=aB_c^b(t) end{aligned}$$
(7)
where parameters a and b are constant in our model and can be used to modulate the plant height growth.
Feedbacks between vegetation and hydromorphodynamics
We consider that the above-ground biomass changes the bed roughness by modifying the Strickler coefficient (K_s), which is used to calculate the flow resistance [Eq. (1)], such as
$$begin{aligned} K_s = K_{s,g}+(K_{s,v}-K_{s,g})frac{B_{c}(t)}{B_{c,max }} end{aligned}$$
(8)
where (K_{s,g}) represents the roughness of the bare bed, which depends on the sediment grain size, while (K_{s,v}) ((<K_{s,g})) is the roughness of a completely vegetated bed assumed to vary with species-specific canopy characteristics. (B_{c,max}) is the biomass density that can be achieved when (B=B_{max}), calculated as (B_{c,max}=lambda _cB_{max}). The presence of vegetation is also known to affect the shear stresses acting on the bed surface and responsible for sediment transport. We model the reduction of bottom shear stress by multiplying the total shear stress (mathbf {tau }) by a factor (gamma <1) and compute the sediment flux using the reduced Shields stress, (gamma theta)32. The parameter (gamma) ranges between 0 and 1 and it is chosen according to (gamma =K_{s}(t)/K_{s,g}), with (K_s) evaluated by Eq. (8).
The role of root-enhanced riverbed cohesion on the evolution of a gravel bed river is taken into account increasing the critical Shields parameter when roots are present. This is implemented as32
$$begin{aligned} theta _{cr}=theta _{cr,g}+(theta _{cr,v}-theta _{cr,g})frac{B_{r}(t)}{B_{r,max}} end{aligned}$$
(9)
in which (theta _{cr,g}) and (theta _{cr,v}) ((>theta _{cr,g})) represent the threshold values for incipient sediment motion on bare and vegetated riverbed, respectively. (B_{r,max}), as (B_{c,max}), represents the biomass density that can be achieved when (B=B_{max}), which is (B_{r,max}=lambda _rB_{max}).
Mortality by burial occurs after the end of any flood event if the amount of deposition in a grid cell is higher than a given fraction of the canopy height, i.e. (Delta z_b > H_{bur}), where (H_{bur}=beta _{bur}H) with (beta _{bur}in [0,1]) (see Fig. 2a) . The parameter (beta _{bur}) accounts for the ability of the plant to withstand sediment burial and the reduction of the canopy height due to bending caused by water flow25. The uprooting is modelled by defining a critical below-ground biomass (B_{r,cr}) that has to be excavated by flow erosion until vegetation is uprooted19. This is defined as32
$$begin{aligned} B_{r,cr}= int _{0}^{hat{zeta }_{upr}}b_r(z)dz = beta _{upr} int _{0}^{1}b_r(z)dz end{aligned}$$
(10)
where (hat{zeta }_{upr}=zeta _{upr}/zeta _r) represents the ratio between the uprooting depth (zeta _{upr}) and the rooting depth (zeta _r). (beta _{upr}) is a constant parameter that defines the strength of the root system to withstand erosion32. Uprooting can occur during the flood event at any time. B and (zeta _r) are set to their initial values when burial or uprooting occurs.
Riverbed erosion and aggradation processes alter the proportion of above and below-ground biomass during floods. According to the mass balance adopted, vegetation is able to re-allocate biomass at a rate that depends on the partitioning coefficients, (lambda _i) [Eq. (6)]. We consider that buried part of above-ground biomass can convert to roots, while exposed part of below-ground biomass caused by erosion can transform in above-ground biomass. Canopy height is then re-calculated using Eq. (7) and the rooting depth is adjusted depending on bed level changes. This assumption is justified by the great plasticity observed in riparian species, which are able to easily resprout from buried stems and grow new tissues from exposed roots49. When no morphological changes occur, the proportion of biomass allocated above- and below-ground can be calculated as (B_i=lambda _iB).
Model workflow
The model workflow is shown in Fig. 2b and considers an alternating sequence of floods and low flow periods. During each flood event morphological changes occur as a result of the two-way interaction between riparian vegetation and river morphodynamic processes. Vegetation can be uprooted during the entire flood event and/or can die by burial at the end of the falling limb of the discharge. Low flow periods are comprised between two consecutive flood events and may last from months to years. In this phase the riverbed is inactive ((theta < theta _{cr})) and vegetation is allowed to grow and develop above- and below-ground traits.
As a first step, we simulate vegetation growth during a low flow period, given a vegetation cover and a riverbed topography. The biomass growth rate [Eq. (4)] is dynamically computed through the evaluation of the function (b_r) [see Eq. (11) in the Appendix], which varies in space depending on the local (cell-wise) variability of the water table level during the growth period. Here we assume that the water table level changes locally with the water surface elevation. This assumption holds true for gravel, uniform substrates where hydraulic conductivity is high43. Discharge variability during low flow periods is responsible for changes in water surface elevations, and thus in water table levels. To evaluate the water surface elevation associated with a certain discharge, we derive a water level-discharge h-Q relation for each computational cell by running one hydrodynamic (fixed bed) simulation for a series of discharges. When the water surface elevation is zero, namely when the cell becomes dry at a certain discharge, we assume that the water table level can be calculated by interpolating the water surface elevations of the nearest wet cells. This allow us to transform discharges in water table level time series and to obtain the frequency distribution of water table levels for each cell during low flow periods. The h-Q relation is then used to calibrate the parameters (lambda _w), (eta _w), (gamma _w) required to compute the function ({b}_r) [see Eqs. (12) and (13) in the Appendix for details] and calculate the growth rate, (sigma _B). By numerically integrating Eqs. (3) and (5) for the duration of the specific growth period, we update the vegetation variables and the associated feedbacks. This procedure introduces a link between plant growth rate and low flow regime characteristics, which in turn influences plant traits.
As a second step, we simulate the riverbed evolution during floods, where feedback mechanisms between vegetation and hydromorphological processes are activated. In particular, the uprooting mechanism [Eq. (10)], the correction of the bed shear stress and flow resistance [Eq. (8)], and the critical Shield parameter [Eq. (9)] are updated every time step. After the flood, we used the new river bed topography for updating the vegetation cover including mortality by burial. The resulting vegetation cover and plant traits are then used as initial conditions for the subsequent growth period.
Source: Ecology - nature.com
