Plant Roots Steer Resilience to Perturbation of River Floodplains

Freshwater ecosystems along river floodplains host among the greatest biodiversity on Earth and are known to respond to anthropic pressure. For water impounded systems, resilience to changes in the natural flow regime is believed to be bidirectional. Whether such resilience prevents the system from returning to pristine conditions after the flow regime changes reverse is as yet unclear, though widely documented. In this work, we show that temporal irreversibility of river floodplains to recover their status may be explained by the dynamics of riparian water‐tolerant plant roots. Our model is a quantitative tool that will benefit scientists and practitioners in predicting the impact of changing flow regimes on long‐term river floodplain dynamics.

with intense riparian vegetation establishment and encroachment causing river channel narrowing (Allred & Schmidt, 1999;Choi et al., 2005;Gordon & Meentemeyer, 2006;Molnar et al., 2008;Stella et al., 2003). From a dynamical system perspective, such ecogeomorphic transformations occur as a "transient phase" that may last for decades (Petts, 1987), before the riverine ecosystem adjusts to a new dynamic equilibrium (Petts, 1984). According to Petts (1987), the "transient phase" depends on several factors including channel type, mobility of sediment and channel boundaries, biota species adaptation, etc. The degree of reversibility of the transformation processes upon restoring pristine hydrologic conditions is largely unknown (Molnar et al., 2008;Perona, Camporeale, et al., 2009;Tullos et al., 2009).
Ecosystem shifts following perturbation have often been ascribed to catastrophe-like dynamics. A tipping point (i.e., bifurcation) toward new stable equilibria occurs when some key system parameter acting as the system driver reaches a critical value (Scheffer et al., 2001). A key feature of such catastrophic transitions is their hysteretic behavior and irreversibility when the system driver conditions are reversed. It is therefore tempting to draw ideas from the catastrophe theory to explain the effect of river impoundment on freshwater ecosystems. However, May (1977) observed that ecosystem dynamics may possess multiple stable equilibrium points and Zahler and Sussmann (1977) pointed out that irreversibility may not necessarily be a consequence of catastrophic transitions. Our present work expands on this idea.
Ecomorphodynamic system theory has elegantly explained how different fluvial styles can be the result of a triad process involving water, sediment, and vegetation dynamics (Bärenbold et al., 2016;Bertagni et al., 2018;Caponi & Siviglia, 2018). A means by which to unravel information and thus quantify the extent and reversibility of floodplain changes to hydrological perturbations is offered by modeling the response of riparian plants and their root systems to perturbation. The analytical tractability of spatial mathematical models inevitably requires simplification even without explicitly considering the dynamics of root adaptation (Bertagni et al., 2018;Caponi et al., 2019). However, further steps in this direction can be achieved by focusing on one point rather than on distributed spatial resolution. This is sufficient to show how resistance to uprooting responds to a changing flow regime and to what extent the process is reversible.
In this work, we develop a comprehensive model that accounts for the evolution of plant uprooting by flow after impoundment, and describes the (stable) equilibrium states of the floodplain system at a point. The complex dynamics of river floodplain response to perturbation are thus reduced to that of a dynamical system represented by a suitable state variable. In particular, we investigate the existence of novel stable equilibrium states for perturbed riverine corridors and discuss their possible irreversibility. The model is applied to a typical example of dam impoundment, which is common worldwide and is known to lead to intense riparian vegetation encroachment with consequent river narrowing (Molnar et al., 2008;).

Riparian Processes and Model Formulation
Figure indicates how roots of phreatophytic vegetation tend to adapt to water table fluctuations. At high elevations above the phreatic surface, the plant root biomass distribution locates preferentially deep into the soil. Conversely, at lower elevations close to the phreatic surface, the root biomass distribution is shallow and highly developed near the soil surface (Tron et al., 2015). Therefore, a vertical (down) shift in the water table may not necessarily hinder the growth of phreatophytic species, but instead affect their rooting depth and vertical distribution . In turn, the anchorage depth of roots influences the ability of a plant to withstand erosion processes and its survival probability of uprooting by flow (Docker & Hubble, 2008;Pasquale et al., 2014;Simon & Collison, 2002). Here, we combine stochastic and deterministic approaches of riparian vegetation dynamics into a comprehensive and almost entirely analytical framework. Accordingly, we use the probability of plant uprooting by flow, P τ as a proxy variable to represent the statistical state of the floodplain at a given time. Hence, it is implicitly assumed that vegetation mortality is solely caused by flow-induced uprooting; other mechanisms such as plant burial have not been considered because their effect may also favor vegetation survival (Kui & Stella, 2016;Politti et al., 2018). Plant uprooting probability depends on plant elevation with respect to riverbed elevation, the representative mean flow erosion event at the riverbed elevation, and the critical scour depth for the plant (Perona & Crouzy, 2018). The second and third factors depend on the statistical properties of river discharge (and water levels), which obviously differ between pre-and post-impoundment periods. We now proceed toward assembling all processes in the model.

Probability Distribution of Time to Uprooting
Perona and Crouzy (2018) modeled plant uprooting by flow as a result of stochastic erosion dynamics requiring a time interval T to scour the bed to the critical depth leading to plant collapse. They obtained the following analytical expression for the probability density function (pdf) of the elapsed time T to uprooting p τ for flow erosion events of generic shape and plant critical rooting depth: where g t (T) describes the noise affecting the erosion process at time t = T, , and τ is the dummy variable of integration. L e is the scouring depth that determines uprooting, and ( ) t is the erosion rate event corresponding to plant elevation. Values of L e and ( ) t are assessed in Sections 2.4 and 2.5. The following expression for g t is obtained assuming that erosion may be represented by an Ornstein-Uhlenbeck stochastic flow process, in which the flow velocity profile is logarithmic and fluctuations acting on sediment particles follow Einstein's diffusion theory (for more details, see the mathematical derivation of Equation 8 in the supporting information [SI]): where D 50 is the median grain size of the sediment, and u * is the shear velocity.

Water Discharge and Groundwater Level Dynamics
Variability in both the water discharge and groundwater levels is addressed using a Compound Poisson process (CPP) (Ridolfi et al., 2011), comprising white shot noise random positive pulses followed by deterministic decays ( Figure 1b). Hence, the pdf of flow discharge is given by BAU' ET AL.

10.1029/2021GL092388
3 of 11 where Q is the flow rate, Γ[⋅] is the Gamma function (Abramowitz & Stegun, 1948), γ d is the mean amplitude of the pulses, and β d is the product between the mean frequency of the jumps λ d , and the deterministic exponential decay rate τ d (see, also SI). Next, we use normal flow conditions to obtain the corresponding water level at each cross section of interest. Likewise, we assume that the water stage follows a CPP with parameters γ l and β l that are fitted to the empirical pdf of the water level (see also SI) and synchronously drive the dynamics of the phreatic water table in the soil (Tron et al., 2014).

Grain Size Distribution
The bed erosion rate and root profile require input values for D 50 , D 10 , and D 90 , which are the median, the 10th, and the 90th percentiles of the sediment size distribution, respectively. To account for the sediment retention capacity of the dam and the reduction in bed mobility downstream, a shift in sediment size between pre-and post-dam periods was included in the modeling framework (Yang et al., 2014). Thus, we did not explicitly model sediment sorting and bed armoring processes; instead, we empirically modeled sediment size increase in the post-dam period.

Root Profile and Scour Depth
According to Perona and Crouzy (2018), the probability of uprooting depends on scouring depth, L e = L 0 − L c , which is the difference between the effective rooting length and the critical rooting length leading to uprooting. We obtain L e by combining the model proposed by Tron et al. (2014) for the vertical root profile r(z), where z is distance below the riverbed level, with that by Bau et al. (2019) for the critical rooting length, c L . As shown in the SI, L e can be obtained by solving the following integral: where L e,t indicates the flow-exposed total rooting length due to scour, and a m is a proportionality constant that links L e,t to its corresponding root biomass, here expressed through the integral of r(z). The mathematical derivation of Equation 4 and further details of the parameters L e,t and a m are given in Supporting Information.

Reference Mean Event and Bed Erosion Rate
Estimation of the probability of uprooting requires knowledge of the temporal evolution of a reference mean erosion event above a given threshold. For the sake of simplicity, we assume that the threshold ξ for onset erosion coincides with the discharge that just starts to inundate the plant at its elevation, η v (see Figure 1). Thereby, erosion (and therefore potential uprooting) at a given location can only occur for flood events whose stage reaches or exceeds the bed elevation at that location, i.e. for values that lie above ξ. To determine the reference flow event, we therefore use the mean of all such events obtained analytically from Calvani et al. (2019), where Q 0 is the mean of all peak events exceeding the threshold ξ, and τ 1 is the integral temporal scale of the reference mean event. The reference mean event (red line in Figure 1c) ceases at   T , which is the up-crossing period of the signal.
From the reference mean discharge event, we then obtain the reference bed erosion event associated with the reference flow event as per Calvani et al. (2019): where λ g is the sediment porosity, ΔX is the erosion length scale, α BL is the coefficient in the bed-load transport formula, A is the wet cross-sectional area of the river, K s is the Strickler coefficient of the sediment, g is the acceleration due to gravity, ρ g is the density of the sediment,  * cr is the critical Shields parameter, and b is the exponent in the sediment transport formula. Equation 6 applies to a point in a generic river section and has been obtained by combining the 1D Exner equation for conditions of net bed erosion (e.g., negligible sediment inflow at the point) with a Meyer-Peter and Müller type sediment transport relationship. The mean erosion event is depicted by the blue line in Figure 1c.

Results From Model Application to an Actual Case Study
The model is applied to the case study of the river Maggia, as it flows through the Valle Maggia in Tessin, Switzerland. After impoundment by dams commenced in 1953, river discharge experienced a severe hydrologic shift which triggered vegetation encroachment and gradual channel narrowing Molnar et al., 2008;. The SI provides a description of relevant data and the calculation of all model parameters. Note that dam impoundment led to a decrease in τ d (from 3.31 to 1.60 days) and λ d (from 0.22 to 0.05 days −1 ), and to an increase in γ d (from 23 to 50 m 3 /s). The product γ d λ d τ d gives the mean flow discharge of the CPP signal μ d . The values obtained for μ d coincide with the mean values of the actual hydrographs, which are 16.5 m 3 /s and 4 m 3 /s for the periods 1933-1953 and 1954-2007, respectively. A table listing the values assigned to the parameters in the equations presented in Section 2 has been included in the SI. Apart from the data retrieved from the literature and previous studies of the Maggia Valley, values of several parameters (related to plant properties and geometry) had to be estimated owing to lack of information. The model satisfactorily represents the expected behavior of the hydrograph, as shown in Figure 2b. During the post-dam period, only the highest river discharge peaks characterize the hydrograph, unlike for the pre-dam period. These peaks correspond to a CPP having higher intensity, lower frequency, and lower temporal correlation. The probability of uprooting P τ was calculated by numerically integrating Equation 1 for the duration of the erosion event and plotting the result as a function of increasing ΔH (i.e., the difference between plant riverbed elevation and mean water stage for increasing hydrograph (down shifts). Hence, ΔH represents the driver in terms of hypothetical hydrologic shift severity caused by impoundment. Figure 2a shows the location of the stable statistical equilibrium states (blue points) of the river floodplain state (represented by P τ ) for increasing ΔH. Pre-dam conditions are represented by the point P1, indicating that plants had more than 90% probability of being uprooted by the reference mean erosion event of the predam hydrograph. The color-rendered aerial photographs (1933,1944, color legend in Figure caption) show the floodplain morphology before dams started to operate. Dam operation produced a vertical downshift in mean discharge (and water stages; see Figure 2b), i.e. a sudden increase of ΔH brought the system to the "out-of-equilibrium" point P2, where plants still had the root architecture of pristine conditions, but were suddenly exposed to post-dam erosion event scenarios. In this case, the probability of uprooting remained high and the image from 1962 shows a floodplain almost without water but with a high braiding index. In the post-dam period (point P2), plant roots started to adapt to the lower water table conditions by deepening root biomass and consequentially reducing the probability of uprooting. This process was gradual, and it took several years for the floodplain system to reach the point P3 (see images 1995 and 2006), which represents the new stable equilibrium for the post-dam hydrological conditions. The same reasoning can be repeated for hypothetical milder shifts of the driver, that is, ΔH, thus obtaining the blue sequence of stable equilibrium states joining points P1 and P3. The process of disconnection between the floodplain and the main channel due to channel incision was not considered here and might lead to two opposite scenarios. For vegetation species able to track the lowering of the water table, their root biomass may deepen further in the soil and the probability of plant uprooting by flooding events would further decrease. On the contrary, vegetation species with low adaptation capability would possibly die and slowly lead to a non-vegetated system, which is not the object of this study. Hence, the first scenario necessarily implies that the value of P τ in the post-dam period may change when considering the ability of different plant species to adapt to extreme and sudden drought conditions. This indirectly explains plant speciation and invasion by species that tolerate and/or favor the new conditions.   1933, 1944, 1962, 1995, and 2006. (b1) Time series of the driver (flow discharge). Note that the value represented for the mean flow rate μ d , is offset for illustrative purposes. (b2) Time series of the ecosystem state (braiding index). The years illustrated are intended to recall the evolution of the spatial distribution shown in the subplots in Figure (a). P1, the deep root system prevents recovery of the original probability of uprooting, thus explaining the tendency of the floodplain to maintain its current narrow morphology.
Given that the model describes only the stable equilibrium points of the system, it is nevertheless instructive to consider the expected dynamics throughout the time domain (Figure 2b2). Up to time t* the system state is at point P1. The state then jumps from P2 to P3 at t = t*, following the hydrologic shift of the driver. From point P2 onward, the probability of uprooting declines, presenting a temporal picture as to how the system states transition from state P2 to P3. The time lapse over which the curve decreases represents the "relaxation time" of the system (in other words, the time required by the ecosystem state to adapt to the new equilibrium). A sensitivity analysis concerning the most relevant input parameters is enclosed in the SI. An important result of the sensitivity analysis emerges when the grain size distribution is maintained constant between pre-and post-dam periods. This preserves the retention capacity of the soil and hence the zone favorable for root growth. For a plant elevation equal to 1.2 m, this results in a value of uprooting probability at the stable equilibrium point P3 that is four times higher than that in Figure 2a. Maintaining sediment continuity in the post-dam period would thus help vegetation control.

Discussions, Implications, and Conclusions
The proposed model has shed light on the type of transitions and temporal irreversibility that potentially affect a river floodplain following hydrological regime shifts. Figure 3 summarizes this process. The curve between P1 and P3 represents the statistical (stable) equilibrium points at which the probability of uprooting follows the progressive adjustment of the root system to the imposed hydrological conditions. In other words, this curve represents a succession of steady states resulting from a quasi-static change in hydrologic conditions. Hence, this segment of the curve can be compared to quasi-static transformations occurring in thermodynamics, where the system always remains at equilibrium. The inability of the system to recover its pristine conditions, for example, such as returning from point P 4 to point P 1 may be ascribed to the development of deep roots as they track the changing water table conditions. However, the model does not explain the dynamical origin of such irreversibility, causing P 4 to also be a stable equilibrium point. This picture appears plausible for water tolerant plants such as riparian plants. Their roots may tolerate long periods under the soil in saturated conditions. Hence, in returning to the original natural flow regime, new deep roots would no longer form, and existing roots might not die off but instead persist in the soil for the entire life time of the plant. Conversely, plant species not tolerating submersed conditions would simply die off and be replaced by others, thus delaying the return to pristine conditions (temporary irreversibility) for which P 4 would be out of equilibrium.
Similar dynamics have also been documented in the literature but have never been modeled quantitatively. For instance, Auble et al. (2007) found that vegetation recovery following the removal of a dam is complex and does not follow a reversal response, leading to the necessity for river restoration intervention. Hence, if dams were to be removed, vegetation coverage and community would not be much affected, leading to a long-term impact on vegetation succession, especially in systems with low sediment transport (Hobbs et al., 2009). The removal of invasive species that often colonize terraces and benches of dammed rivers is extremely complex (Foley et al., 2017), making the process of reintroduction of native species difficult to achieve (Orr & Stanley, 2006;Tullos et al., 2016). For instance, the vegetation response following the dam removal on the Souhegan River in Merrimack (USA) merely consisted of changes to certain herbaceous plants growing closest to the river channel and in the off-channel wetland (Lisius et al., 2018). Furthermore, intensive establishment of mature vegetation during the post-dam period would increase riverbank stability, thus also making it difficult for the river morphology to reestablish its natural pattern (Shafroth et al., 2002). This was also documented by Pearson et al. (2011), who stated that the process of morphological recovery of the Souhegan River has been influenced by the segmentation of alluvial and non-alluvial sections that had been marked by establishment of vegetation on the channel banks during impoundment. Again this, in some measure, is satisfactorily explained by the reduced probability of uprooting by the flow caused by plant root hydrotropic response. Shafroth et al. (2002) also suggested that the persistent occurrence of transient phases after dam construction has a definite impact on the life duration of mature vegetation (e.g., forest), which could persist for even more than a century. In practice, mature vegetation cannot easily be removed by flow erosion processes and return to point P1 may only happen for erosion events of very large return periods or by mechanical action (e.g., restoration). At this point we speculate that a reasonable model representing such out-of-equilibrium system dynamics could have the form where f 1 (P τ ) represents the positive tendency of the system to reduce the root biomass, which would facilitate uprooting. Conversely, f 2 (1 − P τ ), represents the tendency of the system to modify and increase the root biomass in order to decrease the uprooting probability, thus favoring plant survival. Clearly, as P τ depends on the parameter ΔH, the (likely nonlinear) form of f 1 and f 2 should be such that the equation f 1 (P τ ) − f 2 (1 − P τ ) = 0 describes all the equilibrium points (stable and unstable). The fact that the stable equilibrium points of our model joining P 1 and P 3 all lie on a continuous curve suggests that non-reversibility may be ascribed to the presence of other stable equilibrium points (e.g., P 4 ) for the general ecosystem dynamics (May, 1977) rather than to a catastrophic-like mechanism. Such multiple points would represent the capability of water tolerant plants to develop and maintain alive deep roots that tolerate anoxia when conditions are reversed.
Our sensitivity analysis (see the SI) has also shown that the effective particle size of fine sediment plays an important role in the uprooting probability. Hence, replenishment of fine sediment could offer a potential way of maintaining the uprooting percentages for post-dam conditions at levels closer to those for pre-dam conditions. Such a goal could be achieved, for instance, by inducing artificial floods, a well-established technique used to reduce river morphological changes after dam impoundment. In the present application, artificial flooding should be controlled to ensure that the increase in frequency of peak events would bound the erosion rate so as to hinder river narrowing and incision, and stream-bank erosion (Stähly et al., 2019). This strategy could also be adopted to reduce the accumulation of fine sediment upstream of a dam, the presence of which considerably limits the storage capacity of the associated reservoir. The input of fine sediment would not only benefit the shape of the river but also its biodiversity, thus preventing the riparian system from drifting to alternative states (Arheimer et al., 2018). The artificial flooding strategy appears to be promising in terms of effectiveness. This is also confirmed by results obtained by Perona, Camporeale, et al. (2009), who used a lumped model to predict that adding an artificial disturbance each year would lead to increases of 10% in both sediment and water area in the Maggia River reach considered herein.
To conclude, plant root profile can affect riparian ecosystem resilience to pressures such as hydrological alterations and flow erosion processes. Our results suggest that initial state conditions may only be restored after impoundment through the occurrence of a hydrologic event of a much larger return period or by the clearance of riparian vegetation through deforestation and river restoration. This novel combined method can identify and complement dam regulation strategies and promote sustainable solutions to preserve terrestrial and aquatic ecosystems before planetary boundaries are reached (Steffen et al., 2015).

Data Availability Statement
No new data have been used for this research study. Part of the data used to implement the model are available in the following in-text data citation references: Ruf (2007). The aerial images illustrated in Figure 2a can be found online on Mendeley Data (http://dx.doi.org/10.17632/czjwgb9jp8.1).