Bank Erosion Processes in Regulated Navigable Rivers

Vessel‐induced waves affect the morphology and ecology of banks and shorelines around the world. In rivers used as waterways, ship passages contribute to the erosion of unprotected banks, but their short‐ and long‐term impacts remain unclear. This work investigates the effects of navigation on bank erosion along a reach of the regulated Meuse River with recently renaturalized banks. We apply UAV‐SfM photogrammetry, RTK‐GPS, acoustic Doppler velocimetry, aerial and terrestrial photography, soil tests, and multibeam echosounding to analyze the progression of bank retreat after riprap removal. After having analyzed the effects of ship‐generated waves and currents, floods, and vegetation dynamics, a process‐based model is proposed to estimate the long‐term bank retreat. The results show that a terrace evolves in length and depth across the bank according to local lithology, which we clustered in three types. Floods contribute to upper‐bank erosion‐inducing mass failures, while near‐bank flow appears increasingly ineffective to remove the failed material due to terrace elongation. Vegetation growth at the upper‐bank toe reduces bank failure and delays erosion, but its permanence is limited by terrace stability and efficiency to dissipate waves. The results also indicate that long‐term bank retreat is controlled by deep primary waves acting like bores over the terrace. Understanding the underlying drivers of bank evolution can support process‐based management to optimize the benefits of structural and functional diversity in navigable rivers.


Introduction
Human interference on natural dynamics has globally increased to alarming levels, especially on large water courses (Best, 2018). In particular, ship waves can be an important driver of river bank erosion (Houser, 2010;Liedermann et al., 2014;Nanson et al., 1994;Parnell et al., 2007;Styles & Hartman, 2019;Zaggia et al., 2017), with economic (e.g., Rapaglia et al., 2015) and ecological implications (e.g., Gabel et al., 2012). Traditionally, bank erosion is controlled by means of hard measures, such as riprap, but recent approaches searching to balance technical and ecological requirements (Boeters et al., 1997;Heibaum & Fleischer, 2015) allow for some natural processes. Furthermore, in some cases, river banks are allowed to erode again, as for instance along the Meuse River near Gennep, in the Netherlands (Duró et al., 2019). However, banks are allowed to retreat even if there is no comprehensive understanding of ship-induced riverbank erosion and without knowing the final extend of bank retreat.
The knowledge gap is mainly due to the many and complex interacting factors that are involved in the erosion process, especially when ship waves are present. Banks commonly erode because flow currents steepen and undermine them until collapse, after which the toe is temporarily protected by slump blocks (Thorne & Tovey, 1981). Other factors may also destabilize banks such as seepage (Fox et al., 2010) and water-level changes during flood events (Nardi et al., 2012). Ship waves exert additional loads onto banks that may induce mass failures through impinging loads or gradual undermining. Waves can simultaneously act with other drivers enhancing erosion rates. This was observed, for instance, by Dorava and Moore (1997) with currents in bank embayments during peak flows.
The complexity of factors affecting erosion rates involves (i) waves and currents induced by ships that vary in size, speed, loading, and traveling distances from the bank (e.g., Nakos & Sclavounos, 1990), (ii) spatially varying bank geotechnical characteristics (Pollen-Bankhead & Simon, 2009;Samadi et al., 2009), (iii) entrainment rates of bank material , and (iv) vegetation dynamics on eroding banks (e.g., Bertoldi et al., 2011;Edmaier et al., 2011). It is particularly difficult to isolate the effects of the single factors due to their simultaneous occurrence and mutual interactions (e.g., Maynord et al., 2008). In addition, the episodic nature of ship passings hinders the use of survey techniques to measure the effects of single events (Bauer et al., 2002).
One of the aims of river restoration is to increase habitat suitability for fish, invertebrates, and plants by offering a diverse morphology (Wohl et al., 2015) including shallow areas and varied bank slopes. Eroding riverbanks in waterways show a characteristic terrace (Hagerty et al., 1995;Liedermann et al., 2014), which is formed by the combined action of ship-induced waves and stage regulation to facilitate navigation. Banks are hit at approximately the same level over prolonged periods and, as a result, retreat at different rates below and above this level. For this, hereafter, we distinguish the lower from the upper bank ( Figure 1). The lower bank connects the riverbed and the terrace, which are normally submerged due to stage regulation. The upper bank is usually exposed and links the terrace with the floodplain level.
Even though ship waves contribute to terrace formation, increasing the morphological diversity of the river bank, they also negatively affect the local ecosystem through increased bed shear stresses, sediment resuspension, and mobilization of nutrients and chemicals (Gabel et al., 2017). The plants that are able to grow at the upper-bank toe are beneficial for functional diversity (Wollny et al., 2019). Still, their ability to control bank erosion is uncertain (Coops et al., 1996;Koch et al., 2009). A better understanding of bank evolution and long-term response is thus important to define potential ecological improvements, estimate land loss, and evaluate management strategies. The goals of this work are to characterize the processes that determine the evolution of unprotected banks in navigable regulated rivers, integrate them in a conceptual model, and propose a numerical approach to estimate the maximum extension of bank retreat that can be expected. The following aspects are given special attention: relative role of ship waves and floods and role of vegetation growth on bank erosion.

Methods
We analyze bank processes, bank material, ship waves, water flow, and vegetation growth on the left bank of a recently restored reach of the Meuse River, the Netherlands (visible at right-hand side of Figure 2). This river reach is characterized by rather intense ship traffic, presenting a wide variety of erosion rates at the renaturalized bank. For practical reasons, we differentiate here three parts of the vertical bank profile: upper bank, terrace, and lower bank. We study upper-bank erosion processes focusing on the fastest eroding area located within the biggest embayment ( Figure 2). We analyze and compare bank retreat at locations with dissimilar behavior, particularly considering the distribution of erosion before and after a flood event. We analyze the terrace geometry at eight locations, each one having its own physical characteristics. We characterize ship waves based on field measurements and observations. The results are integrated with the aim to define a conceptual model of bank evolution and develop a numerical model to estimate the final bank retreat for different locations.

Study Site
The study site is a 1-km-long reach of the Meuse River located near the city of Gennep, the Netherlands. This site was chosen because it presents different bank erosion magnitudes, which were found to be related to differences in lithological layers laying across the channel (Duró et al., 2019). The Meuse River has a pluvial discharge regime with low flows during the summer, reaching 40 m 3 /s, and floods during the winter, commonly reaching 1,200 m 3 /s and exceptionally 3,100 m 3 /s. The reach under study used to meander before its canalization between 1940s and 1960s, when water levels were also regulated for commercial navigation by means of weirs. Figure 2 shows a 200-m-long and 25-m-deep embayment on the river left bank (observing in downstream direction) during low flows and the minimum regulated water level, corresponding to 8 m NAP (Dutch reference sea level). Positions along the River are indicated in kilometers from the Dutch-Belgium border.

Wave Measurements
We measured water pressure and flow velocity with a Nortek Vector acoustic Doppler velocimeter (ADV) between 21 June and 13 July 2017 with a frequency of 8 Hz. The ADV head and pressure sensor were placed 0.50 m below the regulated water level near the terrace toe where the water depth was 1.60 m (Figure 2). At that location, waves with typical lengths of 3 m were not affected by bottom friction due to deep water conditions. The data were processed with MATLAB scripts. The pressure was converted into water level assuming hydrostatic pressure distribution over the depth. The water-level series was then processed with a third-order median filter to identify primary and secondary waves. Primary waves were isolated by removing fluctuations with frequencies higher than 1/8 Hz and secondary waves by keeping only those frequencies.
The maximum secondary wave height was considered representative for the erosive potential of each secondary wave train (Nanson et al., 1994), so we identified and focused only on those (Sheremet et al., 2013). Appendix A describes the criteria utilized to identify and measure primary waves and the highest secondary wave from the water-level series. summer. Considering only waves with amplitude larger than 5 cm, between 21 June and 13 July 2017, we recorded a total of 1,224 passing ships producing primary waves (Figure 3d), of which 1,013 also produced secondary waves (Figure 3e, circles). In the same period, another 905 ship passings produced secondary waves but induced primary waves smaller than 5 cm (Figure 3e, triangles). Primary wave heights reached 0.45 m, inducing return currents up to 1.2 m/s near the terrace toe. The highest secondary waves also reached 0.45 m, generating orbital velocities up to 0.60 m/s ( Figure 3e). Recreational boats and commercial vessels produced similar secondary waves, but the former usually did not produce significant primary waves (deeper than 5 cm). An example of typical primary and secondary wave patterns produced by a passing ship is available as supporting information Movie S1.
The results confirm that the generation of primary and secondary waves by vessels are not necessarily correlated, given their different geneses (Söhngen et al., 2008). When both types of waves are produced, their relative timing also varies among vessels. The train of stern-secondary waves may happen during the rising limb of the primary wave, as shown in Figures 3b and A1, or after the main primary wave depression (see, for instance, Movie S1). The former case presents coupled water motions consisting of a longitudinal current with overlapping orbital velocities (Figure 3c). These cases allow secondary waves to break at elevations that are lower than regulated levels, that is, over the terrace and below the upper-bank toe (Figure 4a).
At low flows, the breaking location of secondary waves also depends on terrace elevation at the upper-bank toe. If this level is higher than the regulated water level, then secondary waves dissipate on the terrace without reaching the upper-bank toe (Figure 4b, area on the right). On the contrary, secondary waves hit the upper bank when the toe level is submerged or close to the regulated level ( Figure 4b, center left area with vegetation). Secondary waves entrain sediments mostly during breaking, as indicated by the higher concentrations of suspended solids that we observed near the upper-bank toe. Figure 4d shows secondary waves breaking close to the upper bank and previously entrained sediments covering the terrace. Primary waves exert shear stresses during drawdown through the current induced by the transverse energy gradient (i.e., flow toward the main channel), entraining and transporting sediment in suspension ( Figure 4c). The rising limb of deep primary waves turns asymmetric when propagating in shallow water , as above the terrace, with rear slope 2-3 times steeper than the front slope (e.g., Figure 3b). The uprush propagates toward the upper bank as a bore, starting from the lowest water level reached during the depression, subsequently returning the water to the preceding level.  (Figure 9a), when the bankline migrated to form a more uniform retreat from a break at km 153.930, which was caused by a less erodible layer at the upper-bank toe (Duró et al., 2019).

Upper-Bank Erosion Processes
The relative quantities of material eroded during low-flow and flood periods varied along the reach (Figures 5a and 5b). Erosion during low-flow periods is associated to wave action only, given the very low flow velocities at the bank toe, absence of seepage erosion, and rainfall events only producing few local and isolated failures (Duró et al., 2019). During flood events, the flow currents may produce significant shear stresses to contribute to bank erosion (Duró et al., 2019), in addition to ship wave action and water-level fluctuations. From km 153.860 to km 153.885, more bank erosion occurred during periods with flood events than during low flows. From km 153.920 to km 153.955, higher erosion occurred during low flows than in flood periods. Along the rest of stretches, bank erosion happened similar quantities during high-and low-flow periods.  These sections show that the material deposited at the upper-bank toe was not eroded or transported by either currents or waves (km 153.961-153.969, Figure 5b). Near-bank currents were too low at embayments for producing significant sediment entrainment (estimated below 0.3 m/s at the 2018 flood peak while at the channel axis were 1.5 m/s). Waves did not erode the upper bank likely due to the protection of nearby slump-block deposits (km 153.959 and 153.971, Figures 6f and 6i), implying a lateral connectivity and modulation of erosion phases along banks (e.g., see out-of-phase erosion phases of slump-block deposits and undermining in Figure 8b). To conclude, at embayments, floods redistributed the bank material across sections, reshaping the profile with a mild slope (~1:3).

Differences in Bank Retreat
Top-bank retreat shows different behaviors along the studied reach. First, the variety of erodibilities resulted in diverse erosion rates. For instance, km 153.9 evolved twice as fast as 154.1 (Figure 7a) due to different dominant lithology, that is, sandy loams versus silty loams. Moreover, transitions between lithological layers changed erosion rates. For example, rates decreased at km 153.8 at the start of a less erodible layer (Figure 7b), whereas erosion rates increased at km 154.0 after 2014, when a low-erodible layer was surpassed ( Figure 7c, section at pole location).
Uneven distribution of bank retreat also occurred in the largest embayment where floods produced extensive bankline shifts ( Figure 5). Figure 8a shows the topographic changes in that area between the first two measurements, which included a flood event. At the embayment upstream end (around km 153.8), the bank was  (Figure 7b) prevented toe erosion from wave action during regulated water levels. This enabled a rather stable upper-bank slope and profile, also because water-level fluctuations did not destabilize banks once a mild slope was reached, despite the wave attack at higher levels during floods. Similar conditions apply to the upstream ends of the other embayments (Figure 9a), where layers with relatively low erodibility shaped the banklines.
At the downstream end of embayments, banklines crossed layers of different compositions. Here, low-erodibility layers created sharp angles in the bankline that were gradually smoothed by waves and currents ( Figure 9a). Flow currents downstream of protruding banks (see, e.g., km 153.630 in Figure 9a) detached from the upper bank, generating recirculation zones with low velocities over the terrace. At the largest embayment in particular (Figures 2 and 7c), flow currents eventually reattached to the upper bank during floods, due to its length. The reattached currents flowing over the terrace converged at the bay end generating high velocities (up to 1 m/s at the 2018 flood peak). This increased the erosion and transport

10.1029/2019JF005441
Journal of Geophysical Research: Earth Surface capacity compared to more retreated areas, inducing higher shear stresses over protruding low-erodible layers. Nevertheless, bank erosion rates were relatively low, solely attributed to the presence of more resistant layers given the absence of trees at this location.
Vegetation at the upper-bank toe was also found to affect bank erosion rates. Figure 8b shows an area that did not present upper-bank erosion after the 2017 flood because of the presence of sufficiently grown vegetation (km 153.5, Figure 4b, box in Figure 8b). Here, the upper-bank toe remained undisturbed for a sufficiently long time for vegetation to grow thanks to bank retreat rate reduction in 2013 (Figure 7a), possibly due to the encounter of a less erodible layer, and the wave dissipation over a developed terrace. Pioneer plants were first observed in 2015 (Figures 8b and 11b).
The presence of a gravel armoring layer at the upper-bank toe is another factor that affects erosion rates in the study area. Although erosion rates were controlled by the dominant lithology, the gravel layer at upper-bank toe between km 154.275 and km 154.325 ( Figure 4d) reduced the erosion produced by secondary waves after 2013 ( Figure 10d).

Terrace and Lower-Bank Evolution
The terrace topography after 7 years of protection removal presents diverse lengths and water depth at the toe ( Figure 9b) that generally correlate with the different upper-bank retreat rates ( Figure 9a, 2017 bankline). However, each area presented its own variability of terrace toe water depths and lengths. Figures 9c and 9d show the kernel density function of terrace toe water depths and retreats related to subsurface soil cohesion at eight selected areas. The areas are then grouped into low, middle, and high cohesion:  C1, C2, and C3, respectively. Subsurface cohesion shows a general correlation with the two main parameters defining the terrace geometry, since the lowest values correspond to fast retreating areas with the deepest terrace toes, and vice versa.  Figure 10a) present a mildly sloping terrace and considerable upper-bank erosion, whereas high-cohesion areas (C3, Figure 10c) have steeper terrace slopes and low upper-bank erosion. Mid-erodible areas (C2, Figure 10b) show intermediate slopes and upper-bank erosion. Erosion rates across terraces are higher at higher elevations, that is, near the upper-bank toe, and lower near the terrace toe (Figure 10d), except for km 154.3 with an armor layer. Moreover, low-erodible lithologies (C3) present the highest terrace erosion rates, while high-erodible substrates (C1) have the lowest ones.
Different critical shear stresses for entrainment (τ c ) can explain the variety of erosion rates across terrace types. For instance, high-erodible areas (C1) responded faster to wave-induced shear stresses, falling already below τ c in 7 years at the terrace toe, but not at the upper-bank toe where still sufficient wave energy arrives to exert shear stresses above τ c . On the other hand, C3 areas had slow erosion due to high τ c , which combined with still high wave-induced shear stresses result in current relatively high erosion rates. Despite the differences in erosion rates among clusters, terraces evolve in length and depth with a similar spatial and temporal sequence, deepening first the terrace toe and later the subsequent areas across the terrace as it elongates. Final configurations, once shear stresses fall below τ c along the terrace, thus depend on the strength of each lithology. Other factors, however, seem to influence the terrace development. The presence of the armor layer in Area 7 (Figure 9a) can explain the relative short terrace and shallow toe depth compared to Areas 1 and 8 of C2, resulting with similar geometrical characteristics to those of C3 (Figures 9a-9c). Area 6, with a similar terrace geometry to Area 7, is possibly affected by bushes on the floodplain that delay upper-bank retreat. Moreover, Areas 3 and 5 belonging to C3 are both downstream of long embayments, where currents flowing on the terrace converge (see, e.g., Figure 7c and section 3.3) and produce high-flow velocities and thus locally higher erosion rates during floods. It is thus likely that C3 profiles would have shallower and shorter terraces if belonging to homogenous stretches with a uniform exposure to currents.

Ship Wave and Flood Contributions to Bank Erosion
The analysis of data indicates that the terrace is shaped by the regular action of ship waves hitting the bank at regulated water levels during low flows, whereas flow currents are incapable of entraining sediment once a well-developed terrace is formed. Primary waves shear the terrace during drawdown as a current directed toward the channel, and during their rising limb, as a bore traveling on the terrace toward the bank. During low flows, secondary waves regularly act at regulated levels and less frequently at the depression level of the primary waves. During high discharges, secondary waves do not dissipate over the terrace and hit the upper bank.
During early development stages, when the terrace is relatively short, floods appear to contribute to terrace and upper-bank erosion through current-induced shear stresses . This implies flow currents entraining bank material and disaggregating and transporting slump blocks after bank failure (Osman & Thorne, 1988). Furthermore, the presence of ship waves simultaneously attacking banks at high elevations, with low or negligible previous dissipation, promotes further mass failures and block degradation and removal, as observed by Dorava and Moore (1997).

Journal of Geophysical Research: Earth Surface
Such a process is more effective than either factor acting alone. The transport capacity of near-bank currents could move blocks without the need to degrade them (Parker et al., 2011). This is also intensified by the high transverse slopes of the terrace at initial stages (Baar et al., 2018;Thorne & Tovey, 1981). The higher downstream and downslope transport of slump blocks and bank material explains the lower-bank deposits observed between 2010 and 2011 (Figure 10e), when a large flood occurred soon after riprap removal.
Once the terrace develops so that flow currents on terrace and bank are sufficiently reduced (Shiono & Knight, 1991), the flow becomes incapable of removing slump blocks. At this stage, currents and waves decouple their effects on bank erosion and lose their feedback mechanism. At this point, floods operate through water-level fluctuations, destabilizing steep banks particularly during the rising limb of the flood wave. The latter likely happens at the banks of the study site composed by silty clay loam to loam due to loss of negative pore water pressure, with the consequent loss of apparent cohesion and reduction of effective friction angle (Casagli et al., 1999;Duró et al., 2019;Simon et al., 2000;Thorne & Tovey, 1981).
Failed material accumulates at the toe of the upper bank as slump-block deposit. At the same time, waves not only contribute to destabilize steep banks impinging normal forces (Oumeraci et al., 1993) but also rework the slump-block deposits at varying levels during flow recession. These processes of breaking and moving blocks toward the channel are facilitated by their lighter submerged weight. Once the regulated level is reached again, the blocks that are partially or fully submerged on the terrace are sheared by primary and secondary waves. This gradual process, which occurs with mild terrace slopes, progresses until blocks disintegrate and transform in suspended sediment or wash load (ASCE Task Committee on Hydraulics, Bank Mechanics and Modelling of River Width Adjustments, 1998). This mechanism explains the absence of substantial lower-bank deposits between 2013 and 2016 (Figure 10f), despite the continuous upper-bank erosion that occurred during that period (Figure 9a), together with the plumes observed over the terrace.
At this stage, upper-bank erosion occurs with a longitudinal alternation of mass failures and block deposits (herein called modulated failures). This is associated with an increase in form drag during high flows (Leyland et al., 2015), thought to delay erosion rates (Parker et al., 2011). Yet the actual effect remains an open discussion, since it depends on the simultaneous occurrence of modulated failures and the permanence of block deposits. Our observations show the modulation during low flows when slump blocks are not easily removed (see, e.g., Figure 7c), but it may not be the case during flood events where currents clean the toe relatively fast and bank irregularities seem to smoothen during submerged conditions (Konsoer et al., 2017). This would lead to more random and not modulated failures, increasing instead erosion rates by concentrating shear stresses on isolated block deposits.
In deep embayments, variations in volume of eroded material between low-flow and flood periods (see Figure 5) are related to variations in erodibility of the lithological layer that surfaces at the upper-bank toe. Low-flow periods define the duration of the time when ship waves attack the upper-bank toe. The steepness and stability of the upper bank before a flood event depends on previous toe erosion. If the toe is not significantly eroded over 2 or 3 years, as for instance between km 153.840 and 153.860, upper-bank erosion mainly occurs during floods at the bank top level, decreasing the bank slope by mass failures. In this case, the top-bank-level retreats, but not the toe position.
The development stage at which flow and waves decouple their effects depends on three factors. First is river planform, since the highest velocity is found near outer banks (Dietrich & Smith, 1984;Thorne & Hey, 1979). Here, the flow near the upper bank becomes negligible with a longer terrace and vice versa for inner banks. Second is longitudinal flow detachment at outcropping low-erodible layers or trees. This requires a certain length before the shear layer develops and flow reattaches to the upper bank (Van Prooijen et al., 2005), as observed in the longest embayment of the study reach (Figure 2). Third is water depth at the terrace toe, which affects the momentum exchange between main channel and terrace (Knight & Shiono, 1990).

Vegetation and Biofilm Effects on Bank Erosion
Pioneer vegetation growing at the upper-bank toe has been observed at certain locations along the reach 9 years after bank protection removal. It appears that vegetation can only survive the first years of growth under certain conditions. Three locations in the reach presented young trees in July 2019, all of them having a certain degree of protection from wave action, offered for instance, by the presence of rocks on the terrace, which provide extra wave dissipation (km 153.4), or by a groyne placed upstream to stabilize an outlet (km 154.45). In one case, a less erodible layer (intermediate erodibility: Cluster 2) provided both a well-developed terrace and a dry higher ground at the toe, out of the reach of waves during low flows for a couple of years (Figure 4b, km 153.5).
At this location, the upper-bank toe presents an even more resistant soil. Possibly, the availability of a loamy texture also provided favorable physical conditions for plant growth, irrespectively of the necessary preceding seed dispersal (Gurnell, 2014) that naturally happened on the study site. The spatial complexity of propagule dispersal (e.g., Gurnell et al., 2008) could explain the absence of vegetation on other locations with similar conditions to the described above. Km 153.5 particularly did not present adjacent hard structures, and the height of colonizing plants at the toe surpassed the floodplain level (~3 m) in approximately 3 years (Figures 11b and 11c).
The arrival of waves to the vegetated upper-bank toe depended on terrace elevation. Between January 2017 and December 2018, km 153.5 presented a terrace lowering similar to km 154.1 (Figures 10b and 10d). Wave dissipation was progressively less effective over the terrace, increasing the wave energy arrival at the upper-bank toe. Figure 11d shows waves reaching the upper-bank toe where a layer of moderate to low erodibility lays, which holds a sapling on the center right of the photograph. On the center left, a patch with younger vegetation was outflanked by waves, showing higher resistance compared to unvegetated surroundings, and on the other hand, anticipating its removal, as occurred to other young trees in the same stretch (note trunk of dislodged tree on Figure 4c).
It appears that plants could grow on relatively low-erodible soil in the presence of either a well-developed terrace (>12 m) or external wave dissipation, generating a time window without much wave disturbance against the bank. Shorter terraces (Cluster 3) are subject to shear stresses at the upper-bank toe that are too high for plant colonization, whereas longer ones (Cluster 1) present too high erosion rates (due to low τ c ). Low disturbance is necessary for seedling establishment and root growth in other systems too, as for instance on tidal flats (Balke et al., 2011). Even when these conditions occur, subsequent terrace erosion could eventually remove toe vegetation. This agrees with measurements and modeling of salt marshes, which relate marsh retreat to bed-level dynamics of the adjacent flat (Bouma et al., 2016;Mariotti & Fagherazzi, 2013;Willemsen et al., 2018), subject to wind waves among other factors.
Vegetation at the upper-bank toe delays the entrainment phase of the erosion cycle by reinforcing the soil through the root system (Khanal & Fox, 2017). Its presence also protects the upper bank from failing once the repose angle of the slump material is reached (Figure 8; Simon et al., 2011). However, the duration of these effects depends on the terrace stability. Eroding terraces cannot sustain a positive feedback between bank morphodynamics and riparian vegetation dynamics, typical of fluvial and estuarine environments (D'Alpaos et al., 2016;Gurnell & Petts, 2006). Finally, vegetation persistence also depends on the duration of flood events (Glenz et al., 2006).
After 8 years, the terrace presented extended areas covered by biofilms (Figure 11a), likely as a result of shallow water conditions during spring and summer, when high shear stresses are only intermittently induced by ships and light easily penetrates to the bed (Thom et al., 2015). Biofilms reduce the local bed roughness and thus wave dissipation on the terrace. On the other hand, biofilms increase the critical shear stress for sediment entrainment (Cheng et al., 2018;Fang et al., 2014;van de Lageweg et al., 2018), but when this threshold is passed, clumps detach abruptly removing the membrane cover (Vignaga et al., 2013). The penetration of biofilms into the sedimentary bed results in greater erosion resistance over depth, which is sustained over longer time than with superficial layers (Chen, Zhang, Zhou, et al., 2017;Chen, Zhang, Paterson, et al., 2017). Biofilms then affect the development of the terrace (Fang, Chen, et al., 2017;Kakeh et al., 2016;Piqué et al., 2016).

Conceptualization
The factors governing bank erosion in regulated waterways are different for each evolution stage of the terrace. Initial stages are characterized by a relative high contribution of flow currents to terrace and upper-bank erosion and high wave impact. Intermediate stages do not present a significant influence of flow currents on upper-bank dynamics and waves arrive at the upper-bank toe with significant lower energy than in the previous stage. The final stage of terrace development presents a morphology with negligible 10.1029/2019JF005441 Journal of Geophysical Research: Earth Surface upper-bank retreat and terrace erosion, when exerted shear stresses gradually reduce to the critical values for entrainment (final configuration). Figure 12 schematizes the first two settings.
Initial stages (Figure 12a) are characterized by relatively high shear stresses induced by both waves and currents during floods, which drive the upper-bank erosion cycle and terrace erosion. Slump-block dynamics are relatively fast, facilitated by downhill transport. In particular, blocks are generally transported to the lower bank before they are completely disaggregated, contributing to lower-bank accretion. Vegetation could hardly grow in this context, due to either the relatively fast erosion of the upper-bank toe or the frequent high shear stresses at this area. The presence of trees on the floodplain reduces bank retreat rates, but the contribution of trees depends on substrate erodibility at the terrace and upper-bank toe elevation range (Duró et al., 2019).
Intermediate stages (Figure 12b) present lower shear stresses at the upper-bank toe and thus smaller erosion rates, which are more favorable conditions for vegetation growth. Slump blocks are sheared near this area for longer times, due to milder transverse bed slopes and negligible flow currents, disaggregating until the material is entrained, mainly as suspended load. Blocks only contributes to lower-bank dynamics through their sand fraction, which is rather low in the case study, the upper-bank material being mostly silty loam. Small, but continuous, terrace and upper-bank toe erosion results in slow upper-bank retreat, after vegetation decay or removal. Biofilms are able to grow on areas with low erosion rates, temporally influencing critical shear stresses for soil entrainment, but at the same time reducing wave dissipation on the terrace.
Other factors influencing the temporal evolution of the bank profile in the study area are floodplain heterogeneity and the presence of gravel layers. The encounter of a low-erodible layer can drastically change erosion rates, whereas in the opposite case, faster rates are triggered. In the latter cases, longitudinal processes driven by flow currents contribute to smoothen protruding layers. A gravel layer on the upper-bank toe can resist high shear stresses, especially from secondary waves breaking. This creates an armoring effect that delays the upper-bank retreat.
In general, the stability of the terrace depends also on lower-bank stability, which is governed by the shear stresses produced by flow currents and ship return currents. The latter depend on waterway sailing conditions, so, as long as they remain constant, the lower bank could eventually achieve a stable configuration. However, intense floods might always produce erosion. Therefore, the stability of the upper bank depends on terrace stability, which in turn depends on lower-bank stability.
Since the terrace is formed due to the persistent action of ship waves at a close range of water levels, unregulated rivers may have a different morphological evolution due to a load distribution spread over a wider range of water levels. This may also have implications for the long-term development of the bank, especially if vegetation is able to regenerate and reproduce thanks to a low frequency of loads at different elevations. A larger window of opportunity also holds for rivers with a lower ship traffic or milder waves than in the case study (e.g., with larger river cross sections). Nevertheless, the impact of floods under unregulated conditions may be larger due to the proximity and exposure of the bank to currents in the deeper channel. Waves would also not dissipate on a mild sloping terrace, so the net impact on the bank could be even higher. Further research is needed to assess the role of ship waves, floods, and vegetation dynamics in unregulated rivers.

Long-Term Erosion Prediction
The final extension of bank retreat is determined by a balance between the exerted shear stresses and the resistance of both the terrace and the upper bank, if the lower bank remains stable. At minimum regulated water levels, ship waves dissipate on the terrace, exerting progressively lower shear stresses as the terrace extends over time. During floods, both current and waves act directly on the upper bank, destabilizing steep slopes and easily transporting slump blocks away. As the terrace extends, currents gradually lose transport capacity of slump blocks and sediments at the upper-bank toe. The final configuration would therefore suppose a terrace that is long enough to dissipate waves during regulated flows and to substantially reduce near-bank flow during floods, with an upper bank either with a mild (e.g., 1:3) repose slope, or a steeper one if the bank is colonized by vegetation. The presence of biofilm might affect the final configuration by changing the local soil resistance against erosion and the bed roughness, but quantifications of these effects encompass high uncertainties since they depend on several conditions and vary with the season (Fang, Lai, et al., 2017).
To get an insight on the final configuration of the bank, a model is developed to simulate the bed shear stress distribution induced by ship waves at the most unfavorable conditions, that is, the largest wave heights acting at the lowest levels. The model computes the shear stresses induced by primary and secondary waves during their propagation, breaking and running up over a homogeneous terrace, with initial uniform slope. The model updates the terrace morphology with a Partheniades (1965)-type of formula (Equation 1), so that erosion stops when bed shear stresses τ b (Pa) fall below the critical threshold for sediment entrainment τ c . The erodibility coefficient ε is taken as a calibration coefficient. The critical value τ c (Pa) is estimated after Kimiaghalam et al. (2016) formula based on soil cohesion (Equation 2).
Co is the soil cohesion in kPa and ε is the erodibility coefficient in m/s/Pa. The bed shear stresses are computed considering primary and secondary wave action as described below.

10.1029/2019JF005441
Journal of Geophysical Research: Earth Surface

Primary Wave Drawdown
The model considers a constant water-level gradient during the primary wave drawdown, given by the maximum depression H p (m) and the terrace length L (m) to estimate flow velocities through the momentum balance in nondimensional form: Z b is the bed elevation (m), h the water depth (m), U the average velocity over the depth (m/s), g the gravitational acceleration (m/s 2 ), C the Chézy coefficient (m 1/2 /s), and κ = 1/4 an attenuation coefficient (-). The x axis follows the cross-sectional direction and starts at the terrace toe (see later Figure 13a). Since the steady flow assumption with maximum linear water-level gradient was found too conservative to estimate flow velocities, κ compensates for U overestimation, which was estimated by trial and error considering measured velocities. During early development stages, primary waves reach the upper-bank toe before the maximum depression is reached at the terrace toe. As a consequence, the actual maximum energy gradient is lower than the estimated with H p . These cases consider then an effective primary wave height H p,eff (m) instead of H p , computed as where T UBt is the time the primary wave takes to reach the upper-bank toe (s) and T p /2 is half the primary wave period. Friction losses are estimated with the Chézy formula, assuming that the flow resistance at a given section is the one of uniform flow with the same depth and velocity. We adopted the Chézy coefficient given by White-Colebrook formula (Equation 5) and Nikuradse roughness height K s of 0.02 m. The bed shear stress is computed with Equation 6.

Primary Wave Surge
The primary wave propagates as a bore on the terrace during the rising limb. Bore propagation and dissipation is modeled based on energy balance (Battjes & Janssen, 1978, Equation 7), considering refraction and shoaling from an incident angle of θ = 20°normal to bankline (CIRIA et al., 2007) and Snell's law. The water level during propagation is given by the maximum depression reached during the primary wave recession (H p ), assumed horizontal. Since H p was measured from the regulated level, the stern wave (bore) height is H b = 1.5 × H p , taken as upper limit considering measured values and CIRIA et al. (2007) formulation. After breaking, Lamb's (1932) dissipation model for fully developed bores is used (Equation 10).
The primary wave bore is assumed steady with a constant sawtooth shape, resulting in a variance of the wave surface elevation B 0 of 1/12. The bore period (T b ) is taken as 1/3 of the primary wave period (T p ) due to the wave asymmetry in shallow water that steepens the rear slope . T p is taken equal to 25 s for deep waves, which corresponds to the deepest recorded of 0.45 m. We also consider the breaker criterion for nonlinear waves, so γ = H b /h B = 5/9 (Didenkulova et al., 2006;Zahibo et al., 2008). The mean velocity under the trough is estimated through Hansen's (1990) regression for experimental data (Equation 11), and the simplified velocity profile of Svendsen (1984): where β is the ratio between measured and calculated velocities given by with h B being the water depth at breaking. C w is the wave celerity taken as the velocity of the bore (C b , Hansen & Svendsen, 1987), estimated by Svendsen et al. (1978) as with d c being the water depth at the bore crest and d t the water depth at the trough. The bed shear stress is computed with a Chézy-type formula (Equation 14, Jonsson, 1966), with f b a friction coefficient for oscillatory flows computed with Equation 15, the Colebrook formula (Briganti et al., 2018;O'Donoghue et al., 2016).
where Re = u w h/ν is the Reynolds number for instantaneous depth and velocity, with ν being the water kinematic viscosity (m 2 /s).

Wave Runup
Assuming a steady bore, the energy dissipation becomes unrealistic near the upper-bank toe (Putrevu & Svendsen, 1991), for which the formulation changes into runup for h/h B < 0.15. The maximum runup height R v above the location of the bore collapse is estimated with Equation 16  Journal of Geophysical Research: Earth Surface where Fr b,c = U b,c /(gH b,c ) 1/2 is the bore Froude number at the moment of collapse, being U b,c and H b,c mean velocity and height, respectively. α is an empirical coefficient equal to 0.889 accounting for the conversion efficiency of potential to kinetic energy during collapse. The distribution of maximum shear stresses τ b,max is then computed using Pujara et al. (2015) nondimensional relation with the distance along the runup, based on measurements with different types of breaking waves. The upper-limit linear relation of these measurements reads: with x R (m) distance over the runup horizontal length R h (m).

Secondary Waves
The shoaling, refraction, and dissipation of secondary waves are computed with Equations 7-8, with B 0 = 1/ 8 for sinuous waves, and group celerity C g = n × C w . n is computed with Equation 18 and the wave celerity C w = L s /T s through the dispersion relation of linear waves (Equation 19).
L s is the secondary wave length. The wave period T s is equal to 2.25 s and the incident angle normal to bankline θ = 45°, being both the average values from the measurements. The breaker criterion for secondary waves is γ = H s /h B = 0.8, and dissipation is only considered after breaking. Once h/h B < 0.15, runup shear stresses are computed with Equations 16 and 17.

Morphological Update
The initial topography consists of a 1:2.5 slope for the lower bank (s LB ), in correspondence with the 2010 profile of the case study (Figures 10a-10c), with 2.5 m water depth at the downhill boundary. The terrace has a 1/4 × s LB slope, equal to the computed maximum (see below). It starts at a depth equal to 0.40 m, after rounding the deepest recorded primary wave. Nevertheless, the computed profiles have a very low sensitivity to the initial bank slope. The model estimates development stages based on a simplified approach of single-wave heights (one for primary and one for secondary waves), not representing actual erosion rates. In this way, each time step Δt accounts for the erosion of single events, that is, ship passings. After each computational step, the terrace elevation is updated resulting in a new geometry for the input waves. The model is run until no morphological changes occur.
The implicit assumption of the model is that development stages are driven by the frequency of occurrence of given wave heights. The highest waves occur less frequent than the lowest ones (Figure 3), and thus, the terrace morphology takes longer times to adapt to the former. The choice of small input wave heights is representative of initial development stages, whereas higher waves correspond to later stages. Thus, the highest wave happening at a study reach results in the final morphological stage of the bank profile.
The spatial step ΔX for hydrodynamic computations is 0.02 m. The topography is updated through Equation 1 with a morphological spatial step ΔX m of 2 m, considering the mean τ b over −ΔX m /2 and +ΔX m /2 at each point ( Figure 13a). τ b along the terrace results the envelope of maximum values of those induced by the primary wave drawdown, primary wave bore, and secondary waves. Those updated points are then checked for minimum (S min ) and maximum (S max ) slopes, whose values are 0 and 0.10, respectively ( Figure 13a). S min corrects downhill slopes and avoids negative values considering that irregularities are eventually smoothed by erosion, as observed in the field. S max acts uphill and prevents the development

10.1029/2019JF005441
Journal of Geophysical Research: Earth Surface of a scarp, representing the effect of upper-bank retreat without the actual processes involved. The bed elevations of points with ΔX spacing in between those morphologically updated are linearly interpolated.

Modeled Conditions
The model is tested for characteristic values of τ c for Clusters 1-3, respectively, taken as 8, 12, and 18 Pa (Table C1). Two scenarios are modeled for those three cases. In the first, the model is run to represent the profiles measured on 18 January 2017, which are used to estimate ε. The constant-height input wave for each case is evaluated by trial and error within the range of measured waves (Figures 3d and 3e). The waves selected as representative for 2017 profiles are those best reproducing them in terms of terrace length and depth at the toe. In the second scenario, the highest measured waves are used to estimate the maximum bank retreat for each cluster.
The input waves for every simulation are a primary and a secondary wave of independent height per event, representing a ship passing. The number of events to reach an equilibrium morphology is unknown in advance. Every simulation repeats the input primary and secondary waves until the bank morphology does not change between two consecutive time steps. For instance, the necessary number of events for Clusters 1-3 for 2017 profiles was respectively 223, 93, and 77 and for final profiles 434, 291, and 95 (see the modeled profiles later in Figure 13).
The period of primary waves is kept constant and equal to that of the deepest recorded waves, that is, 25 s. The period of secondary waves are set equal to the measured average, that is, 2.25 s. Given that primary waves were found to produce much higher erosion than secondary waves (see next section), only the former were changed among simulations, while the latter were kept with a constant value of 0.45 m, which is the recorded maximum.

Model Results
Predicted 2017 profiles present terrace lengths of 24.9, 12.2, and 7.4 m and water depths at the toe of 0.83, 0.59, and 0.50 m for C1-C3, respectively ( Figure 13, continuous black lines), falling within the measured ranges (Figures 9b-9d). The results indicate that primary wave bores dominate the terrace evolution and final state because they yield the highest shear stresses throughout. This is supported by the fact that secondary waves from recreational boats concentrate in summer, but without an apparent increase of erosion rate over this period ( Figure 5b). The 2017 profiles reasonably match the measured ones with corresponding primary wave heights of 0.30, 0.21, and 0.15 m, which are within the measured wave heights (Figures 3d and 3e). Since the wave frequency of occurrence diminishes the higher the waves are (Figures 3d and 3e), the model shows that the terrace length of each soil type at a certain time, and consequently for a given wave height, depends on soil erodibility, in this case characterized by τ c . More erodible banks have a faster morphological response to higher waves, whereas low-erodible banks need more time to adjust to the same loads. The latter then adjusts to average loads for the given time used to test the model, that is, 7 years. The difference in depth at terrace toes for C2 and C3 can be ascribed to this single-wave height approach that lacks full wave spectra. Such approach could result in a more accurate temporal development, improving the estimate of ε. The effects of currents, also relevant close to the main channel, are not included in the model. Yet the current model is able to represent well the average behavior and characteristic geometry of each cluster.

Discussion of Model Results
The model shows that for the highest measured loads, that is, H p = 0.45 m and H s = 0.45 m (Figure 3), the final development stages are 50%, 100%, and 150% larger than 2017 retreats respectively for C1-C3. This suggests that the lower the erodibility, the slower the terrace evolution. Furthermore, final profiles show that the lower the bank erodibility, the steeper the resulting slope, in accordance with numerical experiments of Bendoni et al. (2019) on tidal flats with cohesive mud-sand mixtures. Field measurements 50 years after stage regulation in the Kanawha River exhibited diverse terrace slopes (Hagerty et al., 1995, Table 1), covering the predicted ones (1:30, 1:21, and 1:17, Figure 13) and exceptionally reaching 1:43 for fine sandy silt, which is comparable to the model estimate of 1:45 for the lowest τ c = 6.4 Pa (Table C1). Figure 14 shows the terrace lengths for τ c and H p ranging 6-18 Pa and 0.15-0.45 m, respectively, which cover the values of the case study, with H s = 0.45 m assuming constant recreational boat waves and T p = 25 s considering the lowest period of primary waves.
Previous models and tools to estimate potential erosion in waterways are rational yet empirical in nature (Glamore, 2008;Spruyt et al., 2012). The model presented here addresses this challenge with a process-based approach, requiring measured wave heights and soil characteristics. The above results correspond to ε = 0.01 m/s/Pa, which was used as calibration parameter to match 2017 terrace parameters. Measured values of ε cover a wide range from 10 −5 to 10 −2 depending on the erosion mechanism, sediment characteristics, and consolidation (e.g., Jacobs et al., 2011;Mitchener & Torfs, 1996). Moreover, the commonly used JET method to estimate ε for consolidated soils encompasses high uncertainty (Karamigolbaghi et al., 2017), and ε is often used for model calibration (Crosato, 2007;). Yet the sensitivity of the model to ε is considerable, especially regarding the water depth at the toe, which for instance reduces to 65% with ε = 0.001 m/s/Pa, whereas the terrace length reduces to 80% with respect to calibrated values ( Figure 13). Nevertheless, the model is able to represent with single-wave heights a coherent response of the system, considering the measured profiles and long-term projections, despite the nonlinearity of the terrace evolution. Future improvements may include wave statistics, which could increase the accuracy of the temporal developments and ε approximations. Other possible improvements are the addition of factors acting in longitudinal direction, such as flow currents, considering upstream detachment or not, and heterogeneous compositions of different layers, including gravel layers. More advance strategies could incorporate the effects of upper-bank erosion processes and vegetation dynamics.

Conclusions
This work aimed to characterize the processes that drive bank erosion in navigable regulated rivers, integrate the roles of relevant factors, and estimate the final extension of bank retreat. We characterized the waves acting on the study site, showing their dynamics and impact locations. We analyzed upper-bank erosion processes, distinct retreat rates across bank profiles, terrace geometry and its relation with soil lithological properties, and the evolution of terrace and lower bank. Finally, we discussed the relative contributions of ship waves and floods to bank erosion, observations on vegetation dynamics. We developed a conceptual model of bank evolution and a numerical model for long-term bank retreat prediction.
Several processes influence the evolution of riverbanks in regulated waterways. The characteristic terrace produced by ship waves attacking at regulated stages develops in two stages, distinguished by the contribution of floods to erosion. Flow currents during floods initially have significant entrainment and transport capacity, but these effects are reduced by the distance of the upper bank to the main channel. In a second stage, floods simply destabilize banks through water-level fluctuations. During low flows and regulated stages, deep primary waves shear the terrace through transverse currents and bores. The latter also erode the upper-bank toe, together with secondary waves that normally attack at this level.
The terrace erodibility given by its lithological characteristics defines the range of lengths and toe elevations at a given time. Soil types have here been clustered in three categories. For each one, the terrace evolves adjusting the length and water depths across the profile, progressively from the toe to the upper bank. Floodplain heterogeneity across single profiles may cause changes from one bank type to another, leading to mixed situations. Currents and waves propagating along the navigation channel tend to smooth bank line transitions. The efficiency of the terrace to dissipate waves, which depends on its length and elevations, eventually controls the upper-bank stability. In turn, the permanence of the terrace position depends on the stability of the lower bank.
Other factors affecting the terrace development include upper-bank dynamics, which present a spatial modulation of out-of-phase undermining and basal clean-out, whose occurrence or net effect with longitudinal currents is not clear yet. Slump-block dynamics also affect the terrace evolution by interacting with dissipating waves and currents, which particularly affect lower-bank dynamics during initial stages. Grown vegetation temporarily protects the upper bank from failure and toe erosion, but its permanence is subject to terrace stability and efficiency to dissipate waves. Moreover, vegetation needs sufficient wave dissipation and dry areas in order to establish at the upper-bank toe. Furthermore, the presence of gravel layers and trees on the floodplain also affects the terrace and upper-bank erosion rates. The long-term terrace stage or final configuration is controlled by the magnitude of primary waves, inducing shear stresses during their propagation as a bore and at breaking, and by the lithological characteristics of the terrace. The necessary time to reach the final stage depends on the abovementioned factors and is mostly influenced by terrace characteristics and ship traffic, as indicated by measurements and suggested by the model results. It appears that the final stage would asymptotically be reached over time, unless biological factors (e.g., biofilms) or anthropogenic interventions change entrainment thresholds or reduce shear stresses. Biofilm effects on cohesive environments should be further investigated to better assess their impacts on long-term morphology, particularly when penetrating into the subsurface.
The development of a terrace across banks of navigable regulated rivers is driven by the additional loads of ship waves acting at minimum stages. The longer the terrace evolves, the more efficiently it dissipates waves, leading to a "self-healing" mechanism. The bank may take decades to reach a stable configuration, in which wave action becomes incapable of entraining more sediment or uprooting vegetation. The factors, mechanisms, and timescales presented and discussed in the present work can help managers in their search to optimize all river functions in future strategies, for instance, protecting against the action of primary waves once a well-developed terrace is reached. Future research including unregulated rivers and noncohesive banks is highly encouraged.

Ship Wave Characteristics A
The quantification of the period and height of primary and secondary waves from the water-level series was done through the identification of key points. First, the water-level fluctuations with respect to the regulated level (elevation zero in Figure A1a) were filtered. The primary waves were separated from secondary waves Figure A1. Quantification of period and height of (a) primary and (b) maximum secondary waves. Figure B1. Location of study reach and analyzed cross sections.

10.1029/2019JF005441
Journal of Geophysical Research: Earth Surface by keeping only frequencies below 1/8 Hz ( Figure A1a). The secondary waves were isolated by keeping only frequencies higher than 1/8 Hz ( Figure A1b).
Primary waves were first identified through the downward zero-crossing at the threshold of −0.05 m from the minimum regulated level (P1). Then, the minimum water level reached by the subsequent trough was determined (P2). This point was used to compute the primary wave height (or depth). After P2, the following upward zero-crossing was identified (P4). Finally, the first upward zero-crossing before P1 was located. The time duration from P3 to P4 was computed to quantify the primary wave period.
The crests and troughs of secondary waves during and after the main primary wave were determined. Then, the highest crest was identified. In the following step, the minimum trough between the previous and the next trough with respect to the maximum crest was determined. The maximum secondary wave height was computed from the difference between the maximum crest and the minimum trough next to it. The period was quantified by multiplying by 2 the time duration between the two points defining the maximum wave height.
Analysed Cross Sections at Study Reach B Figure B1 indicates the location of the study reach, the cross sections measured with RTK-GPS (blue lines), and the cross sections or bank profiles referenced in Figure 5 and presented in Figure 6. Table C1 presents the characteristics of the soil samples taken from the terrace at 10 cm below its surface.

10.1029/2019JF005441
Journal of Geophysical Research: Earth Surface