Parametric transitions between bare and vegetated states in water-driven patterns

Modeling flow variability using a minimalist approach is a challenging task in river science, because of the nontrivial features of the discharge time series $Q\left(t\right)$ (e.g., long-term correlation, intermittent behavior, etc.). The Compound Poisson Process (CPP) provides a parsimonious and robust strategy (24) by addressing the stochastic equation ${\partial }_{t}Q$ = $\mathcal{F}\left(t\right)-Q/\tau$, in which $t$ is time, $\tau$ is the integral temporal scale (namely the integral of the autocorrelation function of $Q$), $\mathcal{F}$ is a shot noise with mean intervals between two pulses equal to $\tau c_v^2$ and mean value $\overline{Q}{c}_v^2$, and $\overline{Q}$ and $c_v$ are the mean and the variation coefficient of the discharge. A sample realization of CPP is reported in Fig. 2A (blue line). At steady state, the PDF of the CPP is Gamma-distributed (Fig. 2 A, Right) and the upcrossing time $T_{\xi }^{+}$—that is, the average time $Q$ stays above a certain threshold $\xi$—is known in closed form (25)where $E\left(a,z\right)$ is the exponential integral function (26). Eq. 1 will be used later on in the triad interaction.

The second element of the triad interaction—sediment transport—is a threshold process that activates when the flow rate overcomes a lower critical value and becomes morphologically effective when an upper statistically rare flow rate $Q_f$ is overcome (Fig. 2A). This so-called formative discharge, an open issue in geomorphology, resets the riverbed by erasing any previous pattern and uprooting all vegetation. The bare flat state so created is the starting point of our analysis.

In a straight river with nonerodible banks, when $Q>Q_f$, the sediment transport morphodynamically triggers the pattern formation of migrating two-dimensional bedforms, called alternate bars (Fig. 1 D and E), which were theoretically studied through linear and weakly nonlinear analyses (27, 28). Most previous studies neglected the role of flow variability, with the notable exceptions of refs. 29 and 30. Starting from infinitesimal perturbations, the bar amplitude $A\left(t\right)$ asymptotically saturates to a finite value $A_s$ in a nonlinear fashion, essentially following a Stuart–Landau equation (Fig. 2C):For our convenience, we define the saturation time $t_s$, such that $A\left(t_s\right)$ = $0.99A_s$ and $A\left(0\right)$ = $0.01A_s$, and we stipulate that a condition for bar formation is that the formative event must have a duration at least equal to the saturation time—that is, $T_{Q_f}^{+}\ge t_s$. An evaluation of the complex growth rate $\mathrm{\Omega }$ and the first Landau coefficient $\sigma$ is here given through the technique of the Center Manifold Projection (CMP), which, unlike the commonly used multiple scaling theory (28), offers a means to find $A_s$ even away from the critical conditions but close to the neutral curve (31, 32).

Let us introduce the following scaling quantities: the river half-width $B$ for the horizontal lengths, the depth $D_0$ for the vertical length, $U_0$ = $Q_f/2B{D}_0$ for the velocity, $B/U_0$ for the time scale, and $u_s$ = $\sqrt{\left(\mathrm{\Delta }-1\right)gd}$ for the sediment velocity ($\mathrm{\Delta }$, $d$, and $g$ being the relative density, particle size, and gravitational acceleration, respectively). The governing dimensionless numbers result: the aspect ratio $\beta$ = $B/D_0$, the relative roughness $d_s$ = $d/D_0$, the sediment transport capacity $\gamma$ = $d_s{u}_s/\left(1-p\right)U_0$ with $p$ the porosity of sediment, and the Froude number $F$ = $U_0/\sqrt{g{D}_0}$. With the above scaling and constant formative discharge $Q_f$, momentum and mass conservation lead to the celebrated dispersion relation for bar linear instability (27) and to the neutral condition ${\mathrm{\Omega }}_r$ = 0 (subscript $r$ refers to the real part). This condition manifests as a marginal curve in the longitudinal wavenumber–aspect ratio plane, with a minimum at the critical point $\left(k_c,{\beta }_c\right)$ (see Fig. 2D). The linear analysis also provides the bar wavelength (fastest growing mode). Further details are provided in Methods and SI Appendix, Section S2.

The weakly nonlinear analysis through CMP leads to Eq. 2, from which one obtains $A_s$ = $\sqrt{{\mathrm{\Omega }}_r/{\sigma }_r}$ and $t_s\sim 6.6/{\mathrm{\Omega }}_r$. Furthermore, as the definition of volumetric discharge reads $Q_f$ = $2\beta F\sqrt{g}{\left(d/d_s\right)}^{5/2}$, both the saturation time and the formative discharge are functions solely of the morphodynamic parameters ($d_s$, $d$, $\beta$, $F$). This allows us to define just two out of the three CPP parameters $\left(\overline{Q},c_v,\tau \right)$, because the third one is given by the above mentioned condition for bar formation—that is, $T_{Q_f}^{+}\left(\overline{Q},c_v,\tau \right)$ = $t_s$. With the aid of Eq. 1, such condition leads to

After the formative event is extinguished, the flow decreases to ordinary values, reducing the water level and letting the crests of the recently formed bars emerge. At this point, vegetation may grow on the dried areas depending on flow variability. Vegetation is assumed to develop by following a logistic Holling type III equation, wherein the carrying capacity, $K^{*}$, can be assumed as a quadratic function of the water depth, $K^{*}$ = $K^{*}\left(D\right)$ (15, 33). In fact, fluvial tree species usually perform the maximum growth when the water table is at an optimal depth: At lower stages, the roots cannot reach the water, while at higher stages, water logging occurs, thus reducing respiration and gas exchange in the root zone (34). In addition, the submerged sites experience a decline in vegetation because of anoxic conditions and uprooting induced by flow drag, which is proportional to the water velocity squared (16, 35). After assuming a linear relationship between flow drag and biomass removal and defining $\psi$ as the density of biomass per unit area normalized to the maximum carrying capacity, we havewhere ${\nu }_g$ and ${\nu }_d$ are dimensionless growth and decay factors of vegetation and are proxies of the plant biomechanical properties in relation to the specific hydromorphological context and $K$ is the normalized carrying capacity (see SI Appendix, Section S4); $\theta \left(\cdot \right)$ is the Heaviside function (i.e., the decay is active only at submerged sites, when the water depth $D$ is positive). We emphasize that, since the water velocity $\mathbf{U}$ = $\left(U,V\right)$ and depth $D$ account for the presence of the bar, the coefficients of Eq. 4 are space- and time-dependent. During ordinary flows, it is likely that sediment transport vanishes in most of the submerged sites, so that the bars can be approximately regarded as stable, until the occurrence of the next formative event. Moreover, the different time scales for bar saturation, order of days, and vegetation growth, order of years, allow us to separate mathematically the process and treat it analytically.

In this description, vegetation develops as a secondary instability (36) over the finite-amplitude bed topography. Encompassing bar topography in space and flow variability in time, a secondary instability is therefore performed by linearizing Eq. 4 about a new basic state (labeled with tilde), namelyIn Eq. 5, the first term r.h.s. refers to the uniform solution, which is dependent on time through the CPP ($D_u\left(t\right)$ is related to $Q\left(t\right)$ through Manning’s equation). The second term accounts for the presence of the finite amplitude bars (see Methods for more details), in which $s$ and $n$ are the longitudinal and transversal coordinates and ${\widehat{D}}_1$ and ${\mathbf{f}}_1$ are the linear perturbation of the water depth and its transversal structure. Similar forms are considered for the other variables. To make the computation analytically accessible, we substitute the discharge stochastic time series with a statistically equivalent periodic one, constituted by a repetition of the typical mean hydrograph. The latter is obtained by the following compatibility conditions (see red curve in Fig. 2A): (i) The duration of the event equals the average interval between subsequent shots, $T=\tau c_v^2$; (ii) the discharge mean value between formative events is preserved; and (iii) the peak value is such that the coefficient of variation $c_v$ of the stochastic series is preserved. In this manner, the time dependency in Eq. 4 becomes periodic and Floquet’s theorem can be applied. Thus, linearizing Eq. 4 around Eq. 5, the solution for the vegetation is $\psi \sim P\left(t\right)e^{\alpha \left(s,n\right)t}$, where $P\left(t\right)$ is a periodic function and $\alpha \left(s,n\right)\in \mathbb{R}$ is the so-called Floquet exponent, which readswhere $\tilde{K}$ = $K\left(\tilde{D}\right)$. We remind the reader that $\alpha$ is spatially distributed and its sign provides the asymptotic behavior of the secondary instability: Wherever $\alpha >0$ vegetation patterns develop. The value of $\alpha$ is more easily evaluated numerically, as the analytical solution of the integral in Eq. 6 is particularly cumbersome. Although the above theory has been described for river systems, it is quite general, and the same steps could be repeated for other ecogeomorphologic systems driven by unsteady flows.

Once the triad interaction is set up (i.e., hydrodynamics, sediment, and vegetation parameters are provided), the present mathematical framework can evaluate the surface area of alternate bars where plants endure flow variability ($\alpha >0$) between two formative events, resetting the morphology. A graphical example of finite amplitude computation and pattern formation is given in Fig. 3A. Let us define the AVI as the ratio of the vegetated area over the emerged area at the minimum discharge, $Q_{min}$—that is, the region theoretically colonizable by vegetation. Fig. 3B shows the dependence of AVI on flow variability $c_v$, for two typical cases: a sand bed river (orange line) and a gravel bed one (light blue line). These two rivers share the same mean discharge $\overline{Q}$, formative event $13\overline{Q}$, and channel width, but they differ in terms of the sediment and the Froude number $F$ (a proxy of slope variation). This, in turn, affects the CPP correlation $\tau$ through Eq. 3. The two cases exhibit analogies and discrepancies. In fact, both rivers show that higher flow variability (increasing $c_v$ and thus $Q_{max}$) reduces the bar portion colonized by vegetation. In particular, there is a threshold value of flow variability, henceforth referred to as $c_v^{*}$, above which plant growth is completely inhibited (namely, AVI decays to 0). In such conditions, vegetation does not have sufficient time to develop since it is frequently removed by the flow. This $c_v^{*}$ threshold corresponds to a parametric transition from vegetated to bare states. Concerning the discrepancies, the longer saturation time of the sand bed river (10 h) with respect to the gravel bed river (2 h) corresponds—via Eq. 3—to a discharge time series that is more correlated in the former case. This leads to longer submergence periods in the sand river and to plant uprooting for a fatigue stress. In the gravel river, the periods of bar submergence are instead shorter, but the more frequent events are equally able to remove mature vegetation. Another discrepancy concerns the two values of $c_v^{*}$, due to distinct slopes and bar heights. Under the same discharge conditions, the gravel bed river exhibits higher bars and shallower water depths and requires the maximum discharge $Q_{max}$ (or flow variability) to be higher to flood the bar top and remove all of the vegetation.

The value of $c_v^{*}$ is a key quantity, which results from an equilibrium between plant growth and decay. It corresponds to the flow variability that nullifies the integral (6) at the most favored site—that is, the bar top. The exact way to compute $c_v^{*}$, which is a cumbersome formulation, is provided in SI Appendix, Section S5. However, as the uprooting mechanism is much faster than plant growth, it is possible to simplify $c_v^{*}$ evaluation assuming that no vegetational pattern develops when the maximum flow ($Q_{max}$) reaches the bar top. In this way, $c_v^{*}$ can be numerically obtained as the solution of the following equation:

Fig. 3$c$ shows that the above approximation is asymptotically exact for very fast uprooting (i.e., ${\nu }_d/{\nu }_g\gg 1$, typical values are of order $10^4$). With the aid of Eq. 7, the behavior of $c_v^{*}$ can be investigated in the parameter space (Fig. 3D). In the region where bedload transport occurs and alternate bars develop, $c_v^{*}$ generally increases when moving from subcritical sand bed rivers ($F<1$, $d<2$ mm) to supercritical gravel bed ($F>1$, $d>2$ mm). This mechanism allows vegetation spread in mountain rivers with high flow variability.

The computation of AVI is here shown for five actual fluvial cases of widespread interest to the scientific community (a summary of river features and validation results is provided in Table 1). The Isère (France) and the Alpine Rhine (border Austria–Switzerland) share almost identical hydrodynamic and sediment parameters (Fig. 1 $d,e$), but flood events of the latter are much stronger and flow variability is higher (37, 38). This implies slightly lower bar elevation in the Alpine Rhine, with frequent inundation and consequent inhibition of plant growth. In contrast, the slightly more regular discharges of the Isère favor vegetation development. Another test case concerns the Vedder canal in Canada (39). Upstream of the canal, the river is braided with multiple vegetated bars. In the canal itself, the bed is narrower with well-developed alternate bars that are easily submerged by ordinary flows. Thus, no vegetation develops in the canal (see SI Appendix, Fig. S5). Finally, the Arc river (France) has undergone a width reduction due to the construction of a highway in the 1990s (40). This triggered a parametric transition between a fully vegetated condition to a weakly vegetated bed (see SI Appendix, Fig. S6), where plants are periodically uprooted by annual flood events.

As also reported in Fig. 4, the present theory quantitatively captures the underlying physics providing a good matching between theoretical predictions and measurements. A partial exception is the actual vegetation cover of the Arc river, which is not fully captured. This might be due to (i) the very strong anthropic influences upstream of the study area (dams and sediment mining) and (ii) the fact that the Arc river has not yet reached the ecomorphological equilibrium and thus we might expect bars to slowly vegetate in the future.

These encouraging results may be useful to depict changes in biogeomorphological styles induced by man in the Anthropocene and of natural origin since the Paleozoic. For example, the model is able to capture the underlying physics for the river environments tested here and provide insights about the physical conditions that must subsist to allow for the colonization and the spreading of plants. Thus, the recent Devonian plant hypothesis, which suggests how plant root evolution might have contributed to the spreading of vascular plants (6), is supported by this theory. In fact, more robust root systems lead to a decrease in the uprooting coefficient ${\nu }_d$ and to a transition to vegetated states (Fig. 3C). However, further implications concerning how established plants have in turn affected landscape evolution are not contained in the present theory. This requires introducing the feedback of vegetation cover on morphodynamics and would explain why our river bar model underestimates the vegetated bar length (Table 1) with respect to the linear theory. Similarly, other processes not accounted for in this analytical theory concern sediment mining, alteration of sediment supply because of hydropower regulation, different root morphologies, and conditions for seed dispersal.