River Ecomorphodynamic Models Exhibit Features of Nonlinear Dynamics and Chaos

Modeling the nonlinear interactions between flow, sediment, and vegetation is essential for improving our understanding and prediction of river system dynamics. Using simple numerical models, we simulate the key flow‐sediment‐vegetation interaction where the disturbance is intrinsically generated by the presence of vegetation. In this case, biomass growth modifies the flow field, induces bed scour, and thus potentially causes vegetation uprooting when erosion exceeds root depth. Our results show that this nonlinear feedback produces deterministic chaos under a wide range of conditions, with complex aperiodic dynamics generated by a period‐doubling route to chaos. Moreover, our results suggest relatively small values of Lyapunov time, spanning 2–4 growth‐flood cycles, which significantly restrict the predictability of riverbed evolution. Although further spatial and temporal processes may add complexity to the system, these results call for the use of ensemble methods and associated uncertainty estimates in ecomorphodynamic models.


Introduction
In the last two decades, vegetation has been recognized as one of the fundamental components of river systems, controlling river evolution at different spatial and temporal scales (Corenblit et al., 2007;Gurnell, 2014).Water flow, sediment flux, and vegetation interact through nonlinear processes that control river system dynamics and generate specific biogeomorphic patterns (Corenblit et al., 2011;D'Alpaos et al., 2016;Murray et al., 2008;Reinhardt et al., 2010).At the landscape scale, vegetation establishment has been shown to promote a transition from multi-thread, bar-braided morphologies to single-thread, meandering channels (Braudrick et al., 2009;Tal & Paola, 2007).At the reach scale, biogeomorphic feedbacks support the formation and evolution of islands and secondary channels (Gurnell et al., 2001), control the dynamics of river bars (Bertoldi et al., 2014;Caponi et al., 2020;Jourdain et al., 2020), and play a critical role in shaping river deltaic marshes (Nardin & Edmonds, 2014), facilitating the formation of drainage channel networks in tidal systems (e.g., Schwarz et al., 2018;Temmerman et al., 2007).Interestingly, the effects of these biogeomorphic feedbacks have been shown to be drivers of hysteresis effects (Bau' et al., 2021) and abrupt shifts (Scheffer et al., 2001;Wang & Temmerman, 2013;Wang et al., 2016), resulting in high sensitivity to rapid forcing and potentially irreversible changes between a state of extensive vegetation encroachment or a completely bare riverbed (e.g., Bertagni et al., 2018;Bertoldi et al., 2014;Corenblit et al., 2011).
Vegetation dynamics are controlled by the balance between vegetation resistance and hydromorphodynamic disturbance, mainly in the form of uprooting by erosion, as this is the only relevant removal mechanism for most fluvial environments (Bankhead et al., 2017;Bywater-Reyes et al., 2015;Edmaier et al., 2011).Bed and bank erosion can occur regardless of the presence of vegetation, as an external morphological forcing; however, vegetation itself produces changes in the flow field, which can lead to local changes in sediment transport vegetation, mediated by intrinisic biogeomorphic feedbacks, exhibits nonlinear dynamics and chaos • For chaos to be possible, the erosion during a flood must be comparable to the root depth • The time scale on which the riverbed/ vegetation evolution becomes unpredictable is short, ranging between 2 and 4 growth-flood periods Supporting Information: Supporting Information may be found in the online version of this article.capacity, thereby inducing erosion and deposition patterns.This process generates an intrinsic feedback whereby the presence and growth of vegetation increases the disturbance (e.g., Bouma et al., 2007;Follett & Nepf, 2012;Le Bouteiller & Venditti, 2014) and thus potentially induces vegetation uprooting.This mechanism directly links the presence of vegetation to the magnitude of disturbance, inducing a negative feedback loop whereby greater biomass causes greater disturbance and thus potentially a decrease in biomass (see Figure 1).This process is common in many fluvial environments where scour occurs in close proximity to individual plants or vegetation patches.Depending on the configuration, scour may occur upstream or laterally of small vegetated patches (e.g., Kim et al., 2015;Tranmer et al., 2024;Zong & Nepf, 2012) or downstream of larger patches covering the entire channel width (e.g., Diehl et al., 2017;Le Bouteiller & Venditti, 2015).
Here we use simplified models that reproduce the main intrinsic interactions between flow, sediment, and vegetation to explore whether vegetated rivers subject to growth and flood disturbance cycles exhibit chaotic behavior (i.e., aperiodic long-term behavior in a deterministic system and sensitive dependence on initial conditions (Strogatz, 2018)) and how this affects the predictability horizon of river morphology.The occurrence of chaotic behavior could modify our approach to predictive ecomorphodynamic models, their development, and the interpretation of model results, as has been suggested for other geomorphic systems (Phillips, 2003;Salter et al., 2020;Stecca & Hicks, 2022).
Numerical models are key for understanding river responses to environmental changes and guiding river management for risk mitigation and ecosystem restoration (Palmer & Ruhi, 2019;Wohl et al., 2015).Understanding biogeomorphic system dynamics is crucial for enhancing ecomorphodynamic model predictions.

Modeling Approach
We employed two distinct models with varying levels of complexity to investigate the potential for chaos occurrence and its characterization.First, we used an ecomorphodynamic model (one-dimensional) that explicitly incorporates the spatial distribution of erosion processes and vegetation biomass.It describes temporal changes in riverbed elevation and vegetation abundance along the streamwise direction (Caponi & Siviglia, 2018).The model incorporates a minimal set of biogeomorphic feedbacks that effectively capture the complex interactions between flow, sediment, and vegetation dynamics.Then we simplified the 1D model by selectively relaxing certain assumptions and removing the spatial component.Both models operate under the key premise that the dynamics of the system are driven by intrinsic feedback mechanisms linking erosion processes, vegetation abundance and mortality (Figure 1). .Vegetation reduces sediment transport, leading to a greater imbalance between the vegetated and bare areas and thus inducing erosion (c).Erosion increases the likelihood of vegetation uprooting, and when scour reaches root depth, uprooting occurs (d).The overall feedback loop is negative: higher vegetation biomass causes greater sediment flux imbalance and more erosion, ultimately resulting in less vegetation.

Hydromorphodynamics
Hydromorphodynamic processes describing the interactions between sediment transport and the flow field are simulated using a one-dimensional approach.The temporal evolution of the riverbed elevation (z B ) is governed by Exner's sediment continuity equation: where p is the sediment porosity and q B is the unit-width bedload flux in the streamwise direction x.We assume that the flow field adapts instantaneously to the riverbed configuration (Parker, 2004), allowing the water surface elevation to be obtained by numerical integration of the gradually varied flow equation.In this way, we were able to integrate the Exner equation with a much larger time step compared to the fully unsteady model, where the hydrodynamics is governed by the Saint-Venant equations.Importantly, this assumption had no discernible impact on the dynamic behavior of our system, allowing a significant reduction in overall computational time.
The global flow resistance is evaluated using the Manning-Strickler method, where the total shear stress is calculated as τ = ρgu|u|/(ks 2 h 1/3 ), where ρ is the water density, g is the gravitational acceleration, u is the vertically averaged flow velocity, h is the water depth, q = uh is the unit-width water discharge, and ks is the Strickler coefficient.q B is evaluated using the Meyer-Peter and Müller formula (Meyer-Peter & Müller, 1948) as a function of the excess of Shields shear stress θ above a threshold θ cr , where θ = τ/(ρ s ρ)gd s and ρ s and d s are the sediment density and diameter, respectively.The governing equations are solved numerically by discretizing the numerical domain with uniformly distributed cells (see Text S1 in Supporting Information S1).

Vegetation Dynamics and Biogeomorphic Feedbacks
Biomass growth in our model is governed by the logistic equation (Bertagni et al., 2018;Caponi et al., 2020), which can be integrated to obtain the following equation: where B is the vegetation biomass, B 0 is the biomass at the onset of the growth period, K is the maximum standing biomass (carrying capacity), σ is the vegetation growth rate, and t is time (Figure S1 in Supporting Information S1).Biomass increases fluid flow resistance by enhancing local roughness, changing flow patterns, and adding drag (Nepf, 2012).In our model we assume that ks varies linearly with vegetation abundance, changing from ks g for bare riverbed to ks v when biomass reaches carrying capacity (Bertoldi et al., 2014;Caponi & Siviglia, 2018).Vegetation also reduces bottom shear stress, affecting sediment transport (Le Bouteiller & Venditti, 2015;Yager & Schmeeckle, 2013).We account for this reduction by multiplying θ by a factor γ (≤1) (Caponi & Siviglia, 2018;Le Bouteiller & Venditti, 2015) in the Meyer-Peter and Müller formula (Meyer-Peter & Müller, 1948).We then use the reduced Shields stress (γθ) to calculate the sediment flux (q B ), where the stress reduction coefficient is determined as γ = (ks/ks g ) 2 < 1. Vegetation uprooting occurs when riverbed erosion exposes a significant proportion of the roots, reducing the anchoring resistance of the plant (Bywater-Reyes et al., 2015;Edmaier et al., 2011;Pasquale et al., 2012).In our model, it occurs when the scour of the riverbed reaches a threshold value ζ upr and is removed on a single cell basis.When vegetation is removed, we set ks to ks g and B to its minimum value B min > 0, a small value to ensure that vegetation grows when starting from a bare riverbed.The value of ζ upr does not change during simulations.We also explored if the addition of further details, such as root-enhanced riverbed cohesion, by increasing the critical Shields parameter, θ cr and dynamic root depth (i.e., ζ upr (t)) affects the dynamic behavior of the system (modeling details are described in the Texts S2 and S3 in Supporting Information S1).

Numerical Simulations
We applied the model to a straight channel with a constant initial slope, in which vegetation can only establish and grow in a limited inner area and occupies the entire width.We chose this configuration to specifically

Geophysical Research Letters
10.1029/2023GL107951 reproduce an erosion process forced only by the presence of vegetation (Caponi & Siviglia, 2018;Le Bouteiller & Venditti, 2014), similarly to the experimental conditions used by Le Bouteiller and Venditti (2014) and Diehl et al. (2017).This choice allows us to isolate the feedback mechanism between vegetation and erosion and thus capture its effects on system dynamics.Although erosion occurs at the downstream boundary in this configuration, it's important to note that the configuration remains applicable regardless of the erosion's location, be it upstream, downstream, or next to small vegetated patches.This setup broadly represents the dynamics at the boundary between vegetated and bare riverbeds, where the presence of biomass intensifies erosion.The main advantage of the chosen configuration is that it provides complete control over the erosion process.
The model is forced by alternating sequences of constant discharge floods and periods of low flow during which vegetation grows (growth-flood cycles) (Figure S2 in Supporting Information S1).We assume that floods repeat deterministically with a fixed period ΔT F and that the riverbed does not evolve during the low flow periods between two successive floods.Since the duration of vegetation growth (months, years) is much longer than the duration of flood disturbance (hours, days), as is typical for growth-disturbance models (Reluga, 2004), we assume that floods occur rapidly and that factors affecting vegetation growth under normal circumstances are temporarily irrelevant during floods (e.g., Bertagni et al., 2018).During floods, the riverbed morphology within the vegetated area adjusts its slope to accommodate the reduced sediment transport capacity due to the presence of vegetation.This leads to erosion and eventually to uprooting (see, Figures S3A and S3B in Supporting Information S1).The final state of the riverbed after a flood event is influenced by the ratio of disturbance to resistance.
A high ratio results in the complete uprooting of the vegetation (Figure S3C in Supporting Information S1), potentially allowing the riverbed to return to its original configuration, especially if the flood duration is sufficient.Conversely, a lower ratio may result in partial (Figure S3D in Supporting Information S1) or no uprooting of vegetation (Figure S3E in Supporting Information S1).In these latter scenarios, the slope within the vegetated patch increases.During low flow periods, vegetation grows logistically toward its equilibrium value (carrying capacity K) (Figure S1 in Supporting Information S1).Vegetation with the highest growth rate, σ = 1, reaches 99.9% of its maximum growth (i.e., B/K = 0.999) at time t = T v .This time scale is used to normalize time as t* = t/T v .Between two consecutive floods, vegetation that survived to the first flood (B ) grows reaching a higher biomass (B + ) just before the second flood (Figure S2B in Supporting Information S1).
Riverbed-vegetation coevolution is controlled by the balance between hydromorphodynamic disturbance and vegetation resistance (Tal & Paola, 2007).We quantify the magnitude of the erosion process (i.e., a measure of disturbance) by the erosion potential E p (Caponi & Siviglia, 2018), computed as the mean value of erosion after the first growth-flood cycle, within the vegetated area with biomass at carrying capacity.Vegetation resistance is quantified by its root depth ζ upr .Then, we can define the ratio between disturbance and resistance W z B = E p / ζ upr ) as a normalized parameter to measure the relative strength of the two components, as suggested by previous work (Bertoldi et al., 2014;Caponi & Siviglia, 2018;Perona & Crouzy, 2018).Vegetation establishment is favored by small values of W z B .Our numerical model explores various W z B values and vegetation growth rates σ maintaining constant flood intervals (ΔT F ).This approach simulates different vegetation types and flood frequencies.Model parameters, typical for gravel-bed rivers, are detailed in Table S1 in Supporting Information S1.

Zero-Dimensional Discrete-Time Growth-Disturbance Model
We developed a conceptual zero-dimensional (0D) discrete-time model (May 1974;Reluga, 2004) with the aim of deriving the simplest model that represents the fundamental physics occurring at the boundary between vegetated and bare areas (Paola & Leeder, 2011).By considering a control volume approach applied to the entire length of the vegetated area (L v ) (Figure S6A in Supporting Information S1), we obtain a simplified description of the erosion and uprooting processes that occur at the boundary between the vegetated and bare areas.Riverbed erosion (E) results from the imbalance between bedload flux entering q IN B ) and leaving q OUT B ) the control volume in a given interval of time (Δt) and is obtained by applying Equation 1 in discrete form: Geophysical Research Letters 10.1029/2023GL107951 In the control volume, the initial slope S is assumed to be constant over time, and the Shields stress used to calculate the entering and leaving bedload fluxes is calculated assuming a wide and rectangular cross-section and uniform flow conditions.It leads, with these assumptions, to: where s is the sediment relative density.The roughness coefficient in the vegetated area and the bedload flux leaving the control volume are computed with the same equation as in the 1D model.In this case, the sediment balance is computed as , which is always negative as γ < 1.Then, the model follows the same workflow as the 1D model (Figure S2 in Supporting Information S1) and retains the same biogeomorphic feedbacks and logistic law of vegetation growth.
The main challenge in deriving the 0D model is the description of the uprooting mechanism.We assume that the uprooting process depends on the ratio between disturbance and resistance W following the jump condition B = U (W)B + .The uprooting function U (W) characterizes the survival of vegetation during floods and exhibits two limiting behaviors.When erosion is negligible and uprooting does not occur (i.e., W → 0), U → 1. Conversely, when erosion removes all vegetation (i.e., W → ∞), U → 0. To investigate the behavior of the uprooting function between these limiting cases, we performed numerical simulations using the 1D model.First, we used a specific configuration that closely resembled the conditions of the 0D model.This approach lumped biomass into a single value, ignoring the spatial distribution of vegetation within the patch (see Figures S4A and S4C in Supporting Information S1).Next, we run the 1D model with a standard setting to determine the effects of the spatial distribution of vegetation on the uprooting function.Our analysis yielded the following uprooting function (Figures S7A and S7B in Supporting Information S1): (5) β ≥ 1 is a nonlinear shape parameter and larger values indicate that even weak erosion can lead to a significant decrease in biomass.Meanwhile, w m is a threshold parameter, representing the minimum value of W for uprooting to occur.These two parameters were obtained by fitting Equation 5 to the results obtained with the 1D model in its full setting (Figure S7C in Supporting Information S1).Remarkably, results consistently demonstrated that β frequently exceeded unity, with a maximum value of 9, while w m ranged from 0 to 0.4 (Figure S7C in Supporting Information S1).The methodology for obtaining the numerical solution of the 0D model is explained in detail in the accompanying Text S4 in Supporting Information S1.

Nonlinear Dynamics Characteristics
Dynamical systems analysis is used to study the asymptotic behavior of the time series of the variables of the 1D model.Results are presented in terms of average vegetation biomass.The same dynamics are observed for the other two variables, that is, bed elevation and water depth.Simulation results show that our 1D model predicts three possible behaviors: (a) constant biomass equal to carrying capacity (Figure 2a); (b) periodic changes between two (as in Figure 2b) or more (but finite) values; (c) aperiodic chaotic fluctuations bounded by a chaotic attractor (Figure 2c).The first behavior is observed when vegetation grows slowly, erosion is weak, and thus no uprooting occurs.Although the vegetation grows over time, the associated erosion at the boundary between the vegetated area and the bare riverbed never reaches the root depth.Thus, the vegetation is free to grow until it reaches the carrying capacity and remains constant over time.The trajectory of the vegetated area in phase space, defined by the averaged values of biomass (x-axis), bed elevation change (y-axis) and water depth (color scale), is a distinct point when the initial growth-flood cycles are neglected and the focus is on the asymptotic behavior (Right Figure 2a).Periodic dynamics occurs when vegetation biomass grows to a constant value during a single growth-flood cycle and is completely or partially uprooted during the following flood.It can take up to three growth-flood cycles for the vegetation to reach a critical biomass that causes significant erosion, leading to biomass removal.The corresponding attractor in phase space is an invariant loop (Right Figure 2b).When the dynamics is chaotic, the vegetation biomass fluctuates, never repeating exactly the same cycles, although the state of the system in phase space is constrained to a chaotic attractor (Right Figure 2c).Typically, vegetation that develops during low-flow periods is partially removed by subsequent floods.In some cases, however, complete removal of vegetation may occur.
"Deterministic chaos" is defined not only by aperiodic fluctuations, but also by sensitivity to initial conditions (Strogatz, 2018;Turchin, 2013).To test this, we perform three independent model runs for each configuration, starting from nearly identical initial vegetation distributions.The results show that stable or periodic trajectories are insensitive to initial conditions and almost indistinguishable (Figures 2a and 2b).On the contrary, when the solution is aperiodic, the three model runs diverge in time, confirming the chaotic behavior (Figure 2c).Furthermore, the three trajectories were used to compute the maximum Lyapunov exponent (λ), as a measure of the divergence of the infinitesimally close trajectories.We used previously published algorithms (Kantz, 1994;Rosenstein et al., 1993) that directly tests for exponential divergence of close trajectories (see Text S5 in Supporting Information S1).As known from dynamical systems theory, positive values of the Lyapunov exponent are characteristic of a chaotic system with aperiodic fluctuations (Figure S5 in Supporting Information S1).In this case, a relevant parameter is the Lyapunov time, defined as the inverse of the Lyapunov exponent and representing the characteristic timescale at which a chaotic system becomes unpredictable.For stable or periodic configurations with trajectories of initially close points that do not diverge, the Lyapunov exponents are negative; the system is predictable and the Lyapunov time → ∞.
To generate a bifurcation diagram, we plot the asymptotic biomass versus vegetation growth rate; it shows how changes in this parameter affect biomass dynamics (Figure 3a).Within the σ-range [0.0, 0.65], solid bands are caused by chaos, and open areas are the result of periodic cycles.A stable period-3 cycle is found in the range [0.22, 0.26], and a stable period-2 cycle is observed in the largest periodic window in the range [0.34, 0.52].For growth rates above 0.69, biomass grows and eventually reaches a periodic orbit.In the chaotic region we consistently found positive values of the maximum Lyapunov exponent.The corresponding Lyapunov time is always equal to only a few (2-4) growth-flood cycles (Figure 3a, red dots and right vertical axis).
Our model exhibits chaotic behavior only when the flood disturbance (and the associated erosion) is of the same order of magnitude as the vegetation resistance (the root depth).When W z B is much smaller than 1, the resistance dominates and the vegetation always reaches the steady state at the carrying capacity (Figure 3b).In contrast, when W z B is much larger than 1, the disturbance dominates and the vegetation is completely uprooted during the flood (Figure 3c).In both these limiting cases, chaos is prevented.
Root-enhanced riverbed cohesion and dynamic root depth do not change significantly the dynamic properties of the system, including the presence of chaotic behavior (Figures S8 and S9 in Supporting Information S1).2b and 2c).

Route to Chaos and Predictability Horizon
While we find compelling evidence for the existence of chaos, the 1D model does not elucidate the mechanisms for the onset of chaos.The bifurcation diagram in Figure 3a does not clearly show a definitive route to chaos.This lack of clarity is mainly due to the uneven spatial variation in biomass distribution and bed erosion within the vegetated area.Therefore, we used the 0D model and analyzed the resulting bifurcation diagram which reveals a clear period-doubling path to chaos (Figure 4a).
The first and second bifurcations occur at about σ = 0.102 and σ = 0.127, respectively, and the chaotic regime starts at about σ = 0.134.The corresponding maximum Lyapunov exponents are negative when the solution is periodic, zero at the bifurcation points, and positive when the system behaves chaotically (Figure 4b).In this case, the Lyapunov time is about 2 growthflood cycles (λ ≈ 0.5), confirming a very short predictability horizon as computed with the 1D model.Despite the further simplification and removal of the spatial component, the 0D model shows a dynamic behavior similar to that of the 1D model, indicating that the biogeomorphic feedbacks and the intrinsic disturbance mechanism are the key ingredients generating the chaotic behavior.
Our results show that the Lyapunov exponent is approximately independent of the model parameters and remains constant when the ratio of disturbance to resistance W is of order one (Figure 4c).Chaotic behavior is prevented when W tends to zero or infinity (Figure 4c, Inset), resulting in negative Lyapunov exponents.In addition, we examine the variability of the Lyapunov exponent as a function of the unit-width stream power, computed as the product of the unit-width discharge (q) by the river slope (S), in a range typical of gravel-bed rivers (Figure 4d).Again, the values are largely constant, between 0.4 and 0.7, and depend more on the choice of parameters β and w m than on the values of the stream power.However, rivers with small stream power (0.02-0.03 m 2 /s) are characterized by smaller values of the Lyapunov exponents and thus longer predictability horizon.

Discussion
Our study reveals that the intrinsic nonlinear feedbacks between vegetation, flow and sediment transport can generate deterministic chaos, thus fundamentally limiting the predictability of the temporal evolution of the system.The results of the proposed models indicate that when chaos occurs, the maximum Lyapunov time is short and can be quantified by a few (2-4) vegetation growth-flood cycles.In natural rivers, the duration of these growth-flood cycles can vary considerably, depending on the type of vegetation (e.g., aquatic vs. riparian) and on the flood magnitude that is able to uproot the vegetation.In low-energy small rivers where macrophytes are the main engineering species (Gurnell, 2014) the growth-flood cycle is likely to be annual.In larger gravel-bed rivers, the return interval of floods that can significantly affect vegetation has been reported to vary from 1 to 2 years for partial vegetation removal on highly dynamic braided rivers (Surian et al., 2015), to several decades for major vegetation renewal (Belletti et al., 2014).Riparian trees grow rapidly, exerting a strong impact on the flow field already after a few years, when the biomass is closer to the ground.Also the resistance of riparian plants develops in short time, where roots of phreatophytic plants grow to reach the ground water level and can easily reach depths >1 m in 1-2 years (Mahoney & Rood, 1998).
Our analysis suggests that chaotic behavior is more likely to manifest when the magnitude of flood disturbance and associated erosion closely matches the resistance exerted by the roots.Therefore, different river systems may exhibit chaotic behavior on different time scales, also depending on the growth stage of the vegetation and its root

Geophysical Research Letters
10.1029/2023GL107951 network.Moreover, our results suggest that external factors that amplify the erosion process, such as bedform migration and bank erosion, are likely to modify the nonlinear dynamics considered here, where the morphological disturbance is uniquely generated by the presence of vegetation.Increased erosion, if not balanced by vegetation resistance, may prevent the emergence of chaos within the system.This intrinsic mechanism is not limited to rivers and can play an important role in different types of vegetated landscapes, where the interplay between time-dependent disturbances and vegetation resistance is key to shaping their evolution, such as in salt marshes, wetlands, and coastal dunes (Bertagni et al., 2018;Goldstein & Moore, 2016;Marani et al., 2010;Schwarz et al., 2018).
Chaotic behavior has been found in fluvial systems considering only sediment-flow interactions, when flow is redistributed in multiple anabranches, like in braided rivers (Stecca & Hicks, 2022) and in deltas (Salter et al., 2020).However, in less complex river morphologies, such as those with migrating bars, the occurrence of chaos remains to be proven, but it is not excluded (Schielen et al., 1993).Our results demonstrate that adding vegetation may enhance the occurrence of deterministic chaos.On the other hand, additional processes, such as flow stochasticity that has a strong control on vegetation dynamics (Bertagni et al., 2018), may add further complexity to the system and can either enhance or mitigate the chaotic behavior resulting from the interactions of flow, sediment, and vegetation.Our modeling approach was not intended to fully characterize a specific river system, but to show that the main intrinsic feedback between vegetation growth and flow disturbance is likely to generate chaos.

Conclusions
We present two simple ecomorphodynamic models that describe the nonlinear interactions between flow, sediment, and vegetation.The main result we obtained is that deterministic chaos is likely to occur when the characteristic erosion depth during a flood is comparable to the root depth of fully grown vegetation.When this condition is met, the predictability horizon is reduced to about 2-4 flood-growth cycles, that is, a decade or less in a typical temperate river.Despite chaotic systems' inherent unpredictability, it is still possible and relevant to use ecomorphodynamic numerical models to quantify their statistical behavior, for example, by using ensemble methods.These methods, commonly used in weather forecasting, involve the generation of multiple realizations of the model (Bauer et al., 2015), while observational data can reduce uncertainty.
Our identification of chaotic dynamics in biogeomorphic systems warrants further investigation.Confirmation of these dynamics across a wide range of rigorously parameterized models may provide robust validation of our findings (Hastings et al., 1993).In addition, the pursuit of direct field observations, although requiring relatively long time series and facing complexities due to flow stochasticity and changes in sediment supply, is increasingly feasible with advanced data acquisition technologies.These technologies offer promising prospects for capturing detailed phenomena related to morphological erosion and above-ground vegetation dynamics.However, obtaining in-depth insights into root structure and plant resistance to uprooting is essential for a comprehensive understanding of biogeomorphic chaos, highlighting an important area for future research.

Figure 1 .
Figure 1.Schematic representation of the biogeomorphic negative feedback loop.Vegetation grows (a) and increases the roughness, resulting in reduced flow velocity within the vegetated area (b).Vegetation reduces sediment transport, leading to a greater imbalance between the vegetated and bare areas and thus inducing erosion (c).Erosion increases the likelihood of vegetation uprooting, and when scour reaches root depth, uprooting occurs (d).The overall feedback loop is negative: higher vegetation biomass causes greater sediment flux imbalance and more erosion, ultimately resulting in less vegetation.

Figure 2 .
Figure 2. Average vegetation biomass over time (left) and in phase space (right) obtained by the 1D model.Three different behaviors are possible (a) stable (b) periodic, and (c) aperiodic dynamics.The parameters used are (a) σ = 0.2, W z B = 1.0 × 10 3 , (b) σ = 0.8, W z B = 2.9, and (c) σ = 0.2, W z B = 2.9.Results from three independent model runs with nearly identical initial vegetation distributions (B1 = 1.0 × 10 5 , B2 = 1.1 × 10 5 , B3 = 1.2 × 10 5 ) show the sensitivity of the solution to the initial conditions.The Lyapunov coefficient is given for the chaotic case for the runs with initial vegetation distribution B1.

Figure 3 .
Figure 3. Bifurcation diagrams obtained by the 1D model for different values of the ratio between disturbance and resistance W z B .(a) W z B = 2.9, chaotic dynamics is widespread (b) W z B = 1.0 × 10 3 , vegetation grows to its carrying capacity and remains constant over time (c) W z B = 1.0 × 10 6 , vegetation is completely removed during floods, resulting in a bare bed after each flood.Biomass before the flood depends on the growth rate σ, and can be close to 0 for very low growth rates or 1 for σ = 1.Blue arrows in panel (a) mark the conditions reported in (Figures2b and 2c).

Figure 4 .
Figure 4. (a) Bifurcation diagram showing total biomass as a function of the vegetation growth rate for the 0D model: the onset of chaos is achieved by classical period doubling.(b) Corresponding plot of Lyapunov exponent.(c) The maximum Lyapunov exponent as a function of the disturbance/ resistance ratio, and (d) of unit-width stream power.