A Numerical Model of Bank Collapse and River Meandering

Meander migration results from the interaction between inner bank accretion and outer bank erosion/collapse. This interaction has been usually treated as a long‐term average of a sequence of erosion events determined by flow hydrographs. Little attention has been paid to the role that individual bank collapse events play on meander evolution. To fill this gap, we developed a numerical model of river meandering that describes explicitly bank collapse. Results show that as bend curvature increases due to meander migration and elongation, the initially scattered locations of bank collapse events converge toward the channel section where bed shear stress attains a maximum. Simulations illustrate the observed catch‐up behavior between inner and outer banks, driven by intermittent bank collapse events. Moreover, bank collapse is found to speed up short‐term meander migration and, consistent with field observations, meanders turn out to evolve toward a state characterized by constant channel width.

and velocity-based bank erosion, implying that individual bank collapse events are neglected. Both assumptions arise from observations showing that meandering rivers, when averaged over the time-scale typical of meander migration, tend to maintain a constant width Lopez Dubon & Lanzoni, 2019).
Recently, studies have been conducted to address the dynamics of each bank independently, and have focused on the relative importance of outer bank erosion versus inner bank accretion (Asahi et al., 2013;Darby et al., 2002;Eke, Czapiga, et al., 2014;Langendoen et al., 2016;Parker et al., 2011;Zolezzi et al., 2012). Although neglecting bank collapse,  suggested the existence of various regimes of bank interaction: both banks eroding, both banks depositing, bar push (faster migration of the inner bank) and bank pull (faster migration of the outer bank), depending on the initial meander configuration and the soil parameters controlling erosion and deposition rates. A river bend can then switch from one regime to another as it evolves toward an asymptotic state where, eventually, both erosion and deposition are roughly equal. More recently, high-resolution field observations have captured a catch-up behavior between inner and outer banks, whereby the river bends widen and narrow in discrete steps, maintaining a statistically steady-state or "equilibrium" channel half width, B eq (Mason & Mohrig, 2019). This provides inspiration to address the paradox of constant channel width, despite the independent and different migration of inner and outer banks.
As for the assumption of a velocity-based bank erosion, simplified representations of bank retreat, carried out at the laboratory scale, have suggested the importance of intermittent bank collapse events. Consequently, the retreat rate is related not only to the near-bank flow velocity, but is also affected by the ratio between bank height and near-bank water depth (Patsinghasanee et al., 2018;Samadi et al., 2013;Zhao et al., 2020).
Here we set up a morphodynamic model involving bank collapse, to describe the catch-up behavior of meander migration associated with the interrelated dynamics of inner and outer banks. Bank collapse is modeled though either a stress-strain analysis, or employing an empirical function derived from laboratory experiments. In order to shed light on the importance of a realistic representation of bank collapse, we focus on meander evolution until incipient cutoff. The results lead us to a conceptual framework of a catch-up behavior that explains the migration of meanders in relation to bank collapse.

Method
The mathematical set-up governing the planform evolution of meandering channels has been thoroughly described (Frascati & Lanzoni, 2013;Seminara, 2006). Here we briefly introduce the ingredients needed to model bank retreat and accretion.
Bank accretion is assumed to be determined only by sediment deposition, whereas bank retreat is given by the sum of flow-induced sediment erosion and gravity-induced bank collapse. Along the banks, sediment erosion and deposition are assumed to act independently and an excess shear stress formula is applied to each process: Here  E and  D are the rates of bank erosion and deposition (m/s), M e and M d are the erosion and deposition coefficients (m/s), both set to 1 × 10 −6 m/s,  b is the near-bank shear stress (Pa),   , c d are the critical shear stresses for erosion and deposition (Pa), which are set to 3 and 10, respectively. R E and R D are proportionality coefficients introduced to probabilistically constrain channel width variations. Indeed, field observations have suggested that rivers tend to maintain a nearly constant width as they migrate laterally Mason & Mohrig, 2019;Nanson & Hickin, 1983). The associated width fluctuations can be described by a specific probability density function (PDF) (Lopez Dubon & Lanzoni, 2019) embedding the intrinsic features of the investigated river reach and the surrounding floodplain. Following Lopez Dubon and Lanzoni (2019), a probabilistic approach is thus employed to control channel width, and a PDF is introduced to modulate the rate of bank erosion and deposition. Specifically, the parameters R E and R D appearing in relations Equations 1 and 2 are calculated according to the cumulative density function of the selected PDF, that is here assumed Gaussian for the sake of simplicity ( Figure 1b). Bank collapse is modeled by means of either a stress-strain analysis or an empirical function developed from laboratory experiments. Stress-strain analysis employs a set of elasticity equations to describe the relation between soil stress and strain within each bank (Duncan et al., 2014). The Mohr-Coulomb criterion and a critical tensile criterion are used to evaluate the stability of the bank on the basis of soil stress (Gong et al., 2018;Zhao et al., 2019). The bank height, H ) is determined by the sum of near-bank water depth, D b , and an upper bank height, H top , computed with respect to the water surface level (Figure 1a). To minimize lateral boundary effects on soil stress distribution, the bank width is set to 5H b for both banks. The bank region is discretized using unstructured cells, with higher resolution near the bank edge. At each time step, the net sedimentation rate over each bank is calculated as the difference   E D ; the computational meshes of each bank are then evaluated and adjusted according to the updated bank profile, defined at several monitoring points.
Alternatively, bank collapse has been modeled as a continuous process (Zhao et al., 2020), relating as follows the contribution of bank collapse to bank retreat, C bc , to the ratio of bank height to near-bank water depth, ZHAO ET AL. where r u is the undermining rate (m/s) of the bank base and r t is the retreat rate (m/s) of the upper bank border. The ratio r u /r t is related to the type of bank failure and soil properties (e.g., soil cohesion and water content). For toppling failures and bare banks made of sand and/or silt, Zhao et al. (2020) suggested to set r u /r t = 0.35. Given that natural riverbanks are commonly covered by vegetation and comprised of clay and silt, we select a relatively larger value (0.7) of r u /r t (Figure 1c).
The erosion rate induced by the near-bank flow is thus amplified by C bc to account for the average effect of bank collapse events, and the bank retreat rate ( R ) is expressed as: Present simulations have been carried out choosing an initially constant reach-averaged channel half width B havg and letting the river to evolve starting from an initial sinusoidal planform (with Cartesian wavelength L b /B havg = 40) and an amplitude A ensuring either small (A = 2B havg ) or moderate (A = 6B havg ) channel axis curvature at the bend apex. We use the same value for morphodynamic and geotechnical parameters in all runs: mean half width to depth ratio  = 27, dimensionless sediment grain size s d = 0.003, shields parameter for the reference uniform flow  u = 0.06, particle Reynolds number p R = 127, soil cohesion  c = 10 kPa, and critical tensile strength  t = 3.5 kPa. Other geotechnical parameters are the same as Gong et al. (2018).
Although it is possible to account for the sheltering effect due to collapsed bank soil by introducing an armoring coefficient (Motta et al., 2014;Parker et al., 2011), we have decided to reduce the number of parameters, embedding armoring effects in the probabilistic approach used to control the channel width (Lopez Dubon & Lanzoni, 2019). Moreover, by integrating the Exner equation over the entire meander length, we account for the continuous variations of the channel slope as a consequence of changes in channel width (Monegaglia & Tubino, 2019). Figure 1d shows the tendency of the reach-averaged channel half width to attain an equilibrium value B eq , when adopting the above described probabilistic approach. In the remaining, we denote as Case A runs performed without bank collapse (i.e., considering only the retreat due to fluvial erosion), Case B runs carried out modeling continuously bank collapse through Equation 4, and Case C runs taking into account the sequence of single bank collapse events through the stress-strain analysis. The suffixes S and M will be used to denote the small or moderate value of the curvature of the bend apex as a result of the amplitude of the initial sinusoidal configuration. Finally, a constant discharge (51 m 3 /s) has been used in the present simulations. Its value has been chosen such that both outer and inner bank are partly inundated by water flow during the initial stages of the simulation.

Results
As bends evolve from an initial sinuous configuration and bend curvature reaches sufficiently high values, bank collapse tends to localize where the bend experiences the maximum bed shear stress ( Figure 2a). The distribution of the positions of bank collapse in the plane (D b /H b , x/L b ) can be fitted through a Gaussian PDF, characterized by a "wide-and-flat" shape (μ = 0.36, σ = 0.19, continuous blue curve in Figure 2a) for small bend curvatures (initial bend amplitude A = 2B havg ), and a "narrow-and-steep" shape (μ = 0.26, σ = 0.03, continuous black curve) for large bend curvatures (A = 12B havg ). A qualitatively similar trend is evident when changing the ratio D b /H b by means of the upper bank height, H top (dashed blue and black lines in Figure 2a): a decreased ratio D b /H b leads to more frequent bank collapses. This typically occurs when bend curvature is relatively low (e.g., during the early stages of simulations carried out with A = 2B havg ) and the collapse location is scattered along the bend, lagging either upstream or downstream of the position of the maximum bed shear stress. Under these conditions, the departures of the bed shear stress from the value characterizing a straight configuration are quite small and therefore, it is the ratio D b /H b that controls bank stability. Bank collapse due to toppling (Zhao et al., 2020) is expected to occur where D b /H b has lower values, leading to a larger bank retreat. On the other hand, as bends evolve or the initial configuration has already a large enough curvature at the bend apex (e.g., for A = 12B havg ), bank collapse generally tends to occur nearby the section where the bed shear stress is maximum, independently of the ratio D b /H b . This occurs because for relatively sharp bends, the bed shear stress at sections characterized by a small ratio D b / H b but with curvatures smaller than that of the bend apex, is too weak to trigger bank collapse, which is consequently dominated by flow-induced erosion rather than bank stability.
Although the inclusion of bank collapse eventually leads to pre-cutoff meander planforms similar to those usually computed by velocity-based, continuous erosion models, the effect on migration rates is evident and depends on bend curvature and how the adopted model accounts for bank collapse. For small bend curvatures (e.g., in the early stages of runs B S and C S ), bank collapse evaluated by stress-strain analysis tends to enhance meander migration, while a slightly lower migration rate is observed when bank collapse is modeled through the parametrization provided by Equation 3. For large bend curvature (e.g., Cases B M and C M , and later stages of Cases B S and C S ), bank collapse increases meander migration, independently of the adopted model. Figure S1 and Text S1 provides a detailed description of the planforms obtained with different bank retreat descriptions. Also, the evolution of channel width toward an asymptotic state depends on the bank retreat model. Although bank collapse eventually leads to an equilibrium width, its effect on the magnitude of erosion/accretion rate is evident. For example, in the presence of small curvatures and modeling continuous bank erosion/collapse, the regime of bank interaction may switch from both banks eroding (early stages of Case B S ) to bank pull (faster outer bank migration) until the asymptotic constant width is reached (Figure 3a). In contrast, a sudden decrease in bank migration is observed immediately after ZHAO ET AL.  The coordinate x is scaled by the Cartesian meander length, L b . The inset is a representation of bank collapse location and of maximum shear stress distribution along the bend. Simulations have been carried out with an initial bend amplitude A changing from 2B havg to 12B havg (B havg = 25 m), and upper bank height H top equal to either 0.5 or 1 m. The fitted curves are normal probability density functions of the position of bank collapses: blue (μ = 0.36, σ = 0.19) and black (μ = 0.26, σ = 0.03) continuous lines correspond to bend amplitudes of 2B havg and 12B havg , respectively; blue (μ = 0.28, σ = 0.1) and black (μ = 0.32, σ = 0.15) dashed lines correspond to an upper bank height of 1 and 0.5 m, respectively. Circles correspond to an upper bank height of 1.0 m, and squares correspond to an upper bank height of 0.5 m. Stars represent the maximum shear stress. (b and c) Plan-and cross-sectional view of bend migration, induced by the interplay between bank erosion, collapse, and accretion. The blue continuous line (t 0 = 3.5 years) and red dashed line (t 1 = 3.6 years) in plot (b) are at the same time as in Figure 3c, when three bank collapse events occur, indicated by the white dashed line in plot (c). The vertical axis ECB represents Elevation with respect to Channel Bottom in front of the bank. Y bc is the retreat driven by an individual bank collapse event, Y e the retreat driven by bank erosion (occurring at every time step), and Y a the advance driven by bank accretion (also occurring at every time step). each bank collapse evaluated by stress-strain analysis, and the erosion rate curve fluctuates as a result of bank collapse events occurring either upstream or downstream (Figure 3b). Since bank collapse is massive and intermittent (compared to the continuous flow-induced bank erosion, see Figure 2b), the migration of meander bends can be described as a cycle consisting of a sudden retreat at the outer bank, followed by a catch-up behavior characterized by the co-occurrence of a weakened outer bank erosion and an enhanced inner bank accretion (Figure 3c). The above catch-up behavior is not only related to a given localized bank collapse, but is also affected by upstream/downstream collapse events. Specifically, the higher the soil cohesion, the less frequent and more massive are the bank collapses, leading to a longer time-scale and smaller catch-up extent (Figures 3b and 3d).   Figure S1e of supporting information). The bars in plots (b and d) correspond to the timing (indicated by x-axis) and scale (indicated by the dimensionless distance Y bc /B havg ) of the local (bend apex, lower gray bar) and adjacent (upper brown bar) bank collapse events. (c) Temporal distribution of the overall extent of bend apex migration due to erosion Y e and accretion Y a . The arrows in plot (c) correspond to bank collapse events occurring at the bend apex (lower gray arrow) and adjacent (one mesh cell) to the bend apex (upper brown arrow). The various quantities on the ordinate axis have been scaled by the reach-averaged half channel width, B havg . The colors of lines are as follows: red, Case A; green, Case B; blue and orange, Case C. curvature (Sylvester et al., 2019), at least for values of B avg /R up to 0.15. The relatively higher migration rate from the upstream inflection point to the apex point indicates that the planimetric pattern of the simulated meandering river is skewed upstream (see Figure S2a of supporting information). The inclusion of bank collapse through stress-strain analysis affects the curvature-migration relation, as a result of the intermittent nature of collapse events (Case C of Figure 4). When bank collapse occurs, the local curvature and migration rate changes abruptly, accounting for the irregular curve between local curvature and migration rate. As bends evolve, the curvature-migration curve becomes more irregular near the apex point (e.g., at year 5), since bank collapse tends to localize where the bend experiences the maximum bed shear stress (Figure 2a). Model results also suggest that meandering channels follow an evolutionary trend whereby  m R first reaches a maximum as the dimensionless local radius of curvature,  R (= avg / R B ) grows, and then decreases with a concomitant increase of  R , thus yielding a hysteresis loop (see Figure S2 and Text S2 in supporting information).

Discussion
The present model suggests that bank collapse, and consequently, maximum bank retreat rates do not necessarily occur where bed shear stress is maximum (Figure 2a), thus challenging the widely adopted linear relation between bank retreat and excess shear stress or flow velocity (e.g., Ikeda et al., 1981). In addition, since bank stability is related to the ratio D b /H b (Simon et al., 2000;Zhao et al., 2020), which constantly changes along the bend and with bend evolution (due to variation in water depth), representations of bank collapse as a function of bank slope and/or upper bank height are debatable when simulating in detail meander morphodynamics. Bank collapse speeds up the migration rate of meanders, enhancing bend amplification and therefore leading to an increase in channel axis curvature. The latter, in turn, regulates bank collapse location. As a bend evolves from an initial sinuous configuration with small curvature, the initially scattered bank collapse events progressively tend to concentrate nearby the section where the bed shear stress is maximum. Our results also provide evidence that for a small-curvature bend, a simplified representation of bank collapse as a continuous process (Equation 3) is inaccurate in the short-term, since meander migration is essentially a discontinuous catch-up behavior (Mason & Mohrig, 2019;Nanson & Hickin, 1983). Therefore, although long-term modeling of meanders can efficiently simulate bank retreat as a continuous process, the short-term evolution cannot neglect bank collapse.
Our study thus highlights the need to describe bank collapse events when simulating meander migration, in order to obtain a quantitatively correct estimate of the rates of bank retreat. The model incorporates the ZHAO ET AL.   (Sylvester et al., 2019). Faster migration rate within the model can be explained by use of a constant flow discharge in the model. In reality, the channel is only occasionally subject to flood events. An intermittency factor (depending on the hydrological regime of the investigated river) has then to be used (Paola et al., 1992) for ensuring a matching between computed and observed migration rates.
interaction between bank collapse and bend planform with a computationally sustainable effort, at least for short-to medium-term predictions. On the basis of the simulations, we suggest that a typical cycle of catchup behavior describes meander migration driven by bank collapse. During Stage I, a single or a series of bank collapses causes the formation of an initial embayment, whose location and size depend on the bend curvature and the ratio D b /H b . At this stage, cross sections affected by bank collapse events experience an increase in channel width. During Stage II, the increased channel width modulates bank erosion/collapse, promotes deposition, and therefore leads to a catch-up behavior between inner and outer banks. The channel width keeps reducing until the deposition rate at the inner bank is small enough. Since the accretion rate of the inner bank is considerably smaller than the retreat rate of the outer bank, a time lag is expected before channel width reaches its equilibrium value. This leads to an increase in erosion rates of adjacent cross sections, and therefore smooths the collapse-induced embayment and speeds up the migration rate of the whole bend. During Stage III, one evolution cycle is completed, and the meander migrates either upstream or downstream, depending on the phase lag between the bend apex and the peak flow location (Lanzoni & Seminara, 2006). As meanders progressively elongate, an asymptotic state may be approached where bank erosion and deposition occur at nearly equal rates, in accordance with the numerical results obtained by . Catch-up behavior has been described by Nanson and Hickin (1983) and Mason and Mohrig (2019) on the basis of field observations. We suggest that meander migration should be interpreted as the interplay between bank pull (bank collapse at the outer bank) and bar push (bank accretion at the inner bank) in discrete steps (i.e., catch-up concept). Because of intermittent bank collapses (Figure 3), the instantaneous outer bank erosion rate is never equal to the inner bank deposition rate over the short-term scale (i.e., several floods), and a catch-up behavior rather than an asymptotic state is a more realistic representation of natural rivers. The channel banks are in general morphologically inactive for long periods of time. They are subject to relevant changes as a result of floods. In particular, bank collapse commonly occurs during the recession period of flood events causing large retreat at the outer bank (Simon et al., 2000), followed by inner bank accretion during periods of below-average flows. When using a constant formative discharge (as assumed in the present study), the effect of a varying hydrological forcing can be accounted for by introducing an intermittency factor corresponding to the fraction of time the channel is actually experiencing flood conditions (Paola et al., 1992;Parker et al., 1998). In the general scenario considered here, this intermittency factor is set equal to 1. Therefore, the time-scale typical of the catch-up behavior within the model (months) is certainly faster than reported in field observations (years, e.g., Mason & Mohring, 2019).
The proposed model is able to capture the peak in migration rate for , a common feature of meandering rivers, but does not describe accurately the subsequent reduction in  m R as the radius of curvature decreases to small values ( Figure S2 in supporting information). A similar result has also been reported by Eke, Czapiga, et al. (2014), who found a weakened tendency of  m R to decrease after reaching a peak value, thus yielding a hysteresis loop rather than a bell shaped curve. This may be due to the limitations of the mathematical model describing the flow field in the bend. For high bend curvatures, the present linearized flow field model is in fact not valid anymore, and the along bend distribution of the secondary currents may vary significantly with respect to mild bends and does no longer vary in the lower range of  R values, owning to the saturation effect described by Blanckaert (2011) and Ottevanger et al. (2013). Further research is thus needed to unravel the interactions between bank erosion/collapse and bend migration, and to fully elucidate the relation    m R R . Finally, the effect of collapsed bank soil is not explicitly accounted for in this model. Clearly, the sediment deposits (slump blocks) consequent to bank collapse may prevent or promote bank erosion (Hackney et al., 2015;Wood et al., 2001). Since deposits of collapsed sediment may play a key role in determining bank erosion rates, research is needed to incorporate their effect in numerical models.

Summary
The interactions between bank erosion and meander migration are essential. Bend curvature and the associated near-bank hydraulics determine the location of bank collapse, which in turn adds to flow-induced erosion, speeding up migration rates and driving changes in bend curvature. As a bend evolves from a small amplitude sinusoidal configuration, the initially scattered bank collapse events concentrate at the outer bank and converge toward the section where the bed shear stress attains a maximum, implying a transition from a bank-stability-dominated state (controlled by the ratio between near-bank water depth and bank height) to a hydraulic-dominated state (dictated by the near-bank distribution of bed shear stress). The present model illustrates the catch-up behavior observed in nature, driven by intermittent outer bank collapse and continuous inner bank accretion. For the same bend, local migration rate increased with local curvature and the inclusion of bank collapse affects the curvature-migration relation, as a result of the intermittent nature of collapse events. Overall, this research improves our understanding of the interaction between bank collapse and meander migration, and therefore provides new insights on alluvial river evolution.

Data Availability Statement
All data shown in the figures of this research can be found in Zhao et al. (2020). The code developed in this study is available on GitHub website (https://doi.org/10.5281/zenodo.4898932).