A Mechanistic Model for Lateral Erosion of Bedrock Channel Banks by Bedload Particle Impacts

Bedrock incision plays a key role in determining the pace of landscape evolution. Much is known about how bedrock rivers incise vertically, but less is known about lateral erosion. Lateral erosion is widely thought to occur when the bed is alluviated, which prevents vertical erosion and deflects the downstream transport of bedload particles into channel walls. Here we develop a model for lateral erosion by bedload particle impacts. The lateral erosion rate is the product of the volume eroded per particle impact and the impact rate. The volume eroded per particle impact is modeled by tracking the motion of bedload particles from collision with roughness elements to impacts on the wall. The impact rate on the wall is calculated from deflection rates on roughness elements. The numerical model further incorporates the coevolution of wall morphology, shear stress, and erosion rate. The model predicts the undercut wall shape observed in physical experiments. The nondimensional lateral erosion rate is used to explore how lateral erosion varies under different relative sediment supply (ratio of supply to transport capacity) and transport stage conditions. Maximum lateral erosion rates occur at high relative sediment supply rates (~0.7) and moderate transport stages (~10). The competition between lateral and vertical erosion is investigated by coupling the saltation‐abrasion vertical erosion model with our lateral erosion model. The results suggest that vertical erosion dominates under near 75% of supply and transport stage conditions but is outpaced by lateral erosion near the threshold for full bed coverage.


Introduction
Bedrock river incision sets the pace of landscape evolution in unglaciated landscapes (Whipple, 2004;Willett, 1999). Bedrock rivers are laterally constrained by rock banks and have intermittently exposed rock beds that incise vertically (Meshkova et al., 2012;. Bedrock rivers form the lower boundary of hillslopes (Perron et al., 2008) and thus are hard points in the landscape that must be cut through to lower the elevation of the whole landscape (Rennie et al., 2018;Venditti et al., 2019). Incision rates of bedrock rivers are commonly modeled as a function of stream power (Anderson, 1994;Hancock & Anderson, 2002;Seidl & Dietrich, 1992;Tucker & Slingerland, 1994;Willett, 1999) or boundary shear stress parametrized from basin slope-area relations (Howard & Kerby, 1983;Howard, 1994;Moglen & Bras, 1995;Stark, 2006;Tucker & Slingerland, 1996;Whipple & Tucker, 1999;Wobus et al., 2006). These models allow for large-scale predictions of landscape evolution over geologic timescales but mask physical processes responsible for bedrock river incision. This makes the predictions of these models difficult to evaluate because the actual erosional processes may differ in important ways from the model assumptions. Process-based models are needed to investigate the relative role of controlling variables such as rock strength, grain size, roughness, water discharge, and sediment supply and to provide more detailed physical explanations (Beer & Turowski, 2015;Huda & Small, 2014;Nelson & Seminara, 2011;Sklar & Dietrich, 2004, 2006Turowski, 2018;Whipple, 2004;Whipple et al., 2000).
Vertical erosion processes are well known, and several models exist to represent them. Whipple et al. (2000) summarized the processes of vertical incision: abrasion by sediment impacts of bedload or suspended load; plucking from the bed by hydraulic forces; chemical and physical weathering; cavitation; and debris-flow scour. Detailed models of the physics of individual incision processes have been developed to predict bedrock river dynamics, including saltation abrasion model (Sklar & Dietrich, 2004); total-load abrasion model (Lamb et al., 2008); plucking model based on the block topple-sliding mechanism (Lamb et al., 2015;Larsen & Lamb, 2016); bedload abrasion, macroabrasion, and plucking model (Chatanantavet & Parker, 2009); and weathering model (Hancock et al., 2011). These models have been used to predict how vertical incision in bedrock channels changes in response to changing boundary conditions (Egholm et al., 2013;Huda & Small, 2014;Larsen & Lamb, 2016;Sklar & Dietrich, 2006Whipple, 2004).
However, bedrock rivers can also erode laterally and adjust their width. Undercut walls are evidence of active, local width adjustment (Figure 1). Local variations in bedrock river width can induce highly turbulent plunging flow as water enters the narrow part of bedrock rivers, which can in turn promote erosion of the bed and sidewalls by bedload particle impacts (Venditti et al., 2014). Lateral incision has also been observed to be responsible for formation of strath terraces (Fuller et al., 2009), creation of wide valley bottoms (Snyder & Kammer, 2008) and planation of valley bottoms (Cook et al., 2014) at large scales. Therefore, understanding lateral erosion mechanisms is crucial for exploring bedrock width dynamics and its influence on fluvial processes from local (reach) to large scales. In comparison to what is known about vertical erosion, however, comparatively little is known about lateral erosion mechanisms. Previous studies mostly relate the bedrock bank erosion to local conditions, such as flood events (Stark et al., 2010), high alluvial cover (Gilbert, 1877;Shepherd, 1972), meander migration (Finnegan & Dietrich, 2011), weak lithology (Montgomery, 2004;Stark et al., 2010), or bank strength (Limaye & Lamb, 2014). Existing lateral erosion models rely on the stream power law to link stream power or parametrized shear stress to erosion rates with various degrees of sophistication (Croissant et al., 2019;Finnegan et al., 2005;Hancock & Anderson, 2002;Lague, 2010;Langston & Tucker, 2018;Stark, 2006;Turowski et al., 2009;Wobus et al., 2006;Yanites, 2018). Most of these models ignore the influence of sediment supply on lateral erosion by simply scaling the lateral erosion rate with shear stress (e.g., Stark, 2006;Wobus et al., 2006) or the rate of energy dissipation per unit area of the channel wall created by centripetal acceleration around a bend (Langston & Tucker, 2018). Others have introduced the influence of alluvial cover on limiting lateral erosion in high sediment supply environments (Hancock & Anderson, 2002;Lague, 2010;Yanites, 2018) but did not include a quantitative relation between sediment supply and lateral erosion rate because the fundamental relation was unknown. Turowski (2020) recently developed a lateral bank erosion model due to bedload particle impacts, deflected by gravel bars. The model does not include the physics of deflections but rather treats the gravel alternate bars as a source of roughness capable of deflecting particles in an otherwise straight bedrock channel. In the limit of small degrees of cover, this produces decreasing lateral erosion rates with increasing extent of alluvial cover because gravel bars increase their length as the cover gets greater due to the assumption of constant aspect ratio of gravel bars. Gilbert (1877) first suggested that a bedrock channel will incise laterally when the channel bed is covered with transient alluvial deposits. Recent research on lateral erosion has focused on the role of sediment supply on setting the relative rates of vertical and lateral erosion (Finnegan & Balco, 2013;Fuller et al., 2009;Turowski et al., 2007). These investigations suggest lateral erosion dominates in high sediment supply environments but is limited in low sediment supply environments. None of these studies propose a specific process or mechanism to explain how the high sediment supply drives lateral erosion. Physical experiments have documented channel widening by bedload abrasion (Finnegan et al., 2007;Johnson & Whipple, 2010). Enlighted by these experiments, Fuller et al. (2016) further explored the erosional mechanism of deflection of saltating bedload particles into the channel wall by roughness elements and concluded that it is an effective mechanism for lateral erosion into bedrock when the bed is covered by alluvium or composed of large protruded roughness elements. This mechanism explains why lateral erosion dominates in high sediment supply environments where intermittent alluvial cover likely occurs. The downstream transport of bedload particles is deflected by alluvial cover and obtains lateral momentum to erode the wall. In low sediment supply environments, alluvial cover may not be available to deflect bedload particles. This newly identified mechanism for lateral erosion opens the door for a mechanistically based lateral erosion model.
Here we develop a mechanistic model to explore the potential efficacy of bedload particle impacts as a mechanism of lateral erosion in bedrock channels and test the model using the Fuller et al. (2016) flume experiments, referred to as Fuller experiments hereafter. Our model only considers the collision between bedload particles and bed roughness elements as the sole process by which saltating bedload particles obtains lateral momentum to erode the wall. We acknowledge that other lateral erosion mechanisms certainly exist, such as plucking (e.g., Beer et al., 2017), but abrasion is the dominant process in massive crystalline rock. We also recognize that channel curvature may enhance particle impacts with the wall (Cook et al., 2014;Langston & Tucker, 2018;Mishra et al., 2018;Turowski, 2018) but have elected not to include that effect to enable a solution to our model. The numerical model is formulated by determining the initial velocity of bedload particles before collision with bed roughness elements from empirical relations (Sklar & Dietrich, 2004), estimating the momentum transfer during collision from a simplified reflection methodology, and tracking the movement of bedload particles from collision with bed roughness elements to impact on the wall using force balance equations. This allows the distribution of lateral erosion on the wall to be calculated. The model is implemented with and without coevolution of wall morphology, shear stress, and erosion rate to explore how channel change influences the results. After we show how our lateral erosion model works, we then couple it with the Sklar and Dietrich (2004) vertical incision model to investigate the competition between vertical and lateral erosion with transport stage and relative sediment supply.

Model Development
The model is based on the saltation-abrasion mechanism of bedrock erosion and the well-known tools and cover effect (Sklar & Dietrich, 2004). Erosion rates are a function of sediment supply, transport stage, grain size, and rock strength (Sklar & Dietrich, 2004. When the bed is relatively free of cover, impacts of saltating bedload particles are capable of detaching rock particles from the surface. Vertical erosion is limited at high sediment supply rates, when the bed is covered. However, when covered, downstream transport of saltating bedload particles can be deflected by particles that make up the alluvial cover (referred to as roughness elements) and directed toward channel walls, which induces lateral erosion. Following the saltation-abrasion vertical erosion model formulation (Sklar & Dietrich, 2004), we assume that the flow, sediment transport, and distribution of roughness element are uniform in a bedrock channel with a planar bed and straight walls. We use a hybrid approach to model lateral erosion by impacts of saltating bedload particles. First, we model all the possible individual deflection trajectories from discrete parts of the roughness elements for a given hydraulic condition. Then we apply these results in a continuum model by calculating the deflection rates on each cell of the roughness surface and calculate the resultant erosion rates as a function of locations on the wall.

Initial Hydraulic, Flow Resistance, and Bedload Transport Conditions
We assume that bed roughness elements are composed of immobile semispheres with diameter of D r and an areal fraction of F r , arranged in uniformly distributed rows and columns with a spacing of d ( Figure 2). The distance between the wall and the center of the first roughness element is the same as the spacing between two adjacent roughness elements (Figure 2), for simplicity. We tested the effects of this simplification by setting distance to the wall from the first defector to d − D r and found it had a negligible effect on the erosion rates (<1%). Initial hydraulic conditions are calculated from six input variables: water discharge Q w , channel width W, channel slope S, areal fraction F r of roughness elements, roughness element diameter D r , and bedload particle diameter D. In natural bedrock rivers, D and D r would be the grain size of the transported bedload and deposited bed material, respectively.
Assuming steady uniform flow, the total shear stress τ is given as where ρ w is water density, g is gravity acceleration, and h is water depth.
τ can also be expressed as a function of Darcy-Weisbach hydraulic friction factor f and mean flow velocity U: In a bedrock channel with roughness elements and transported bedload particles, the flow resistance is derived from the bedrock surface, roughness elements, alluvial cover, and channels walls. To calculate the contribution of each source, flow resistance has been weighted by its areal proportion (Ferguson et al., 2019;Inoue et al., 2014;Johnson, 2014;Tanaka & Izumi, 2013). Here we adopted the Johnson (2014) method and assumed the wall flow resistance is negligible, which is valid for a channel that is wide relative to its depth. f can be expressed as a weighted average of the spatial fractions of different sources of flow resistance in the channel: where F a is the fraction of alluvium and f b , f r , and f a are friction factors for bedrock, roughness elements, and alluvium, respectively. Because the deposition of alluvial cover was observed to be negligible in the Fuller experiments, equation 3 for that case can be simplified to f b and f r can be modeled using appropriate roughness length scales in any preferred flow resistance relation. For simplification, they are used here as fitting parameters to calibrate the model to the Fuller experiments.
Combining equations 1-4 with the continuity equation for a rectangular channel (Q w = WhU), h, U, and τ can be solved as Assuming the roughness elements cause flow separation and contribute form drag, the shear stress available to transport sediment τ s can be obtained from replacing f in equation 2 with (1 − F r )f b : Figure 2. (a) Cross-section view and (b) plan view of model setup in an idealized rectangular channel eroded by saltating bedload particles that are deflected by roughness elements distributed on the channel bed. The gray semispheres represent roughness elements with diameter of D r , which are equally distributed in rows and columns with the same distance d. The green spheres represent bedload particles that impact roughness elements. Only one side of the channel walls is shown here and used for simulation, assuming that the walls are symmetrical.
Initial bedload transport conditions, including the saltation hop height h s , saltation hop length l s , and bedload particle velocity u s , are estimated from the empirical relations of Sklar and Dietrich (2004): where ρ s is the sediment density, u * ¼ ffiffiffiffiffiffiffi ffi ghS p is the flow shear velocity, τ * s ¼ τ s = ρ s − ρ w ð Þ gD is the nondimensional shear stress available for sediment transport, τ * c is τ * s at the threshold of motion for particle movement, and w f is the particle fall velocity, which is calculated from the empirical method developed by Dietrich (1982), assuming values of Cory shape factor (0.8) and Powers scale (3.5) typical for natural gravel grains. Auel et al. (2017aAuel et al. ( , 2017b recently updated the Sklar and Dietrich (2004) equations but calibrated hop lengths, heights, and velocities to shear stress for supercritical flows with Froude numbers up to 5. The relations notably produce symmetric particle trajectories, which are different from the low velocity, asymmetric trajectories typical of subcritical and transcritical flows we consider here. Application of our model to in supercritical flows may require reparameterization.
The bed-normal velocity w s is calculated from the difference between the gravitational acceleration of the particle and deceleration due to drag (Lamb et al., 2008): where is the drag coefficient, and h c is the height above the bed of the point of collision with the roughness element (h c = 0 for collisions with the bed).

Collision Between Bedload Particles and Roughness Elements
Assuming that the saltating bedload particles have negligible lateral momentum during the normal course of a downstream hop, the saltation lateral velocity v s before collision is zero. Thus, the incoming saltation velocity vector i s has two nonzero components: During collision with roughness elements in water, bedload particles experience an inelastic rebound that can be modeled by the sum of an elastic and a viscous force (Cundall & Strack, 1979). For simplicity, the elastic response is modeled using a reflection methodology to calculate the outgoing saltation velocity vector after collision with a roughness element as where p is the projection of the incoming particle velocity vector onto the surface normal vector, at the point of collision (defined by the normal vector b n) calculated from assuming that the tangential force during collision is negligible. The coefficient of restitution (C r ) describes the retention of particle momentum during the collision between bedload particles and roughness elements. We choose a value C r = 0.9 based on experimental observations (Joseph et al., 2001;Joseph & Hunt, 2004;Niño et al., 1994;Schmeeckle et al., 2001), which means that the particle loses 1 − C 2 r ¼ 19% of its incident kinetic energy during an impact.
The magnitude and direction of o s = (u o , v o , w 0 ) are controlled by incoming velocity i s and normal vector b n at the point of collision (Figure 3). Consider a bedload particle that collides near the base of the roughness element, at the roughness element centerline. The magnitude of i s for this case is maximized because the collision occurs near the bed where w s is the greatest, which will maximize the magnitude of o s for given hydraulic conditions. However, the collision will create an o s for this case that points in the upstream direction with negligible lateral velocity v o , because b n is pointing upstream (Figure 3). In contrast, when b n is rotated to 45 degrees relative to the centerline of the roughness element (Figure 3), o s will have a substantial wall-normal velocity component v o with negligible downstream velocity component u 0 . Therefore, to incorporate the variation of i s , b n, and hence, o s at the point of collision with the roughness element, the surface of each roughness element is discretized into N approximately uniform triangular grid cells (N ≈ 2,000 is selected here for a balance of efficiency and accuracy). Within each cell, the potential impact position and impact angle are assumed to be represented the cell centroid, and the outgoing velocity o s of individual bedload particles is calculated in each grid cell from equations 13-15.
Not all cells on the surface of a semispherical roughness element are subject to collisions. To estimate which cells will experience collisions, and the impact rate on each grid cell as a function of the bedload flux, we begin by assuming that the trajectory of bedload particles before impacting on the roughness element is composed of two components: upward trajectory and downward trajectory (Sklar & Dietrich, 2004). The upward trajectory has a hop height of h s and a hop length of l su , and the downward trajectory has a hop height of h s and a hop length of l sd . Assuming these two trajectories together form a triangle, with a total hop length of l s and hop height of h s (Figure 4), l su and l sd can be approximated from l s as (Sklar & Dietrich, 2004)   Three planes are formed by the triangular trajectory of bedload particles: 1) the plane parallel to the upward trajectory; 2) the plane parallel to the downward trajectory; and 3) the plane parallel to the bed ( Figure 4). All upward-moving particles must move parallel to the first plane, and all downward-moving particles must cross the first plane. In contrast, only upward-moving particles will cross the second plane, and all downward-moving particles will follow the second plane. The third plane, the channel bed, is where the particles turn around. Our model only incorporates the impacts of downward-moving particles on the roughness element surface. The length L of the first plane for intercepting the downward-moving particles is and its angle α intersecting with the bed is The impact rate, with dimensions of collisions per unit time per unit area on the first plane, can be expressed as where q s is sediment supply per unit width and M is the mass of a bedload particle. The area of each grid cell is projected onto the first plane, along a vector parallel to the downward trajectory of bedload particles, to calculate the impact rate on each grid cell. The angle β of the projected direction intersecting with the bed is ( Figure 4) The projected area for each grid cell is defined as A c . The impact rate on each grid cell of the roughness element surface I c can hence be expressed as The variation of impact rates I c is illustrated in Figure 5, where the 10-mm roughness element from the Fuller experiments is used as an example. The center of each grid cell is projected onto a horizontal 2D surface from the plan view ( Figure 5a) and onto a vertical 2D surface from the along-stream view ( Figure 5b). There is no impact on most of the downstream-facing part of the semisphere surface (Figures 5a-b) because it is below the tangent point of downward-moving trajectory. Meanwhile, the impact rate is zero near the vertex of the upstream-facing part of the roughness element (Figures 5a-b), because the impacts here are in the shadow of downward-moving trajectory when the roughness elements are too close to each other ( Figure 4). The impacts decrease from the center to the edge of the roughness element ( Figure 5a), due to the decrease of the shadow effect as the radius r of a circle for a longitudinal slice through the sphere reduces to zero at the edge of the roughness. The impact rate also decreases with distance downstream because the impact area A c goes to zero when the surface cell gets tangential (parallel) to the flux trajectory ( Figure 5b).

Movement of Bedload Particles From Collision With Roughness Element to Impact on the Wall
After collision, the movement of bedload particles is modeled from force balance equations and tracked over each time step Δt. We assume that fluid drag and gravity are the dominant forces affecting instantaneous downstream velocity u, lateral velocity v, and vertical velocity w. The changes in particle velocities with time are given by where U z is the downstream flow velocity at height z above the bed. For turbulent boundary layer flow in a channel, U z can be calculated from the law of the wall where κ is von Karman's constant (~0.41), k s is the hydraulic roughness length scale which can be obtained from friction factor f using a general Manning-Strickler formula (Johnson, 2014) Equations 23-25 can be numerically integrated at each time step Δt to solve for the velocity and position of individual bedload particles. The time step used in the simulation is Δt = 10 −5 s. Smaller time steps were also tested, which substantially increase the computational time but do not change the results. A minimum wall-normal velocity v min is adopted here to distinguish between impacts that cause erosion and impacts that are viscously damped, which is a function of the particle Stokes number S t (Davis et al., 1986;Joseph & Hunt, 2004;Schmeeckle et al., 2001): where ν is the kinematic viscosity of the fluid (10 −6 m 2 s −1 ) and a value of S t = 100 is selected here from Schmeeckle et al. (2001) and Joseph and Hunt (2004). At each time step, a bedload particle may be rebounded by the channel bed or other roughness elements before it impacts on the wall ( Figure 6). In this situation, the rebounded velocity is simulated using the same method used for the original collision with the roughness element, taking into account that the normal vector for the bed is vertical. The simulation runs until a bedload particle has impacted the wall or its lateral velocity is below the velocity limit v min before reaching the wall. When bedload particles impact the wall, the impact velocity vector IV = (u i , v i , w i ) and impact position vector IL = (x l , y l , z l ) are recorded for calculation of lateral erosion rate of different locations on the wall.
The deflection trajectories of bedload particles vary with the impact positions on the same roughness element. Bedload particles impacting on the part that is near 45 degrees relative to the centerline of roughness element travel a shorter downstream distance because the particles have larger lateral velocity and can impact on the wall faster ( Figure 6a). Meanwhile, bedload particles deflected by the higher part of the roughness element can impact higher on the wall due to the higher initial height before deflection and the upward-moving velocity after deflection here ( Figure 6b). When the roughness elements are located further from the wall, more impacts are viscously damped and are rebounded by the bed before impacting on the wall due to more loss of momentum on the way to the wall ( Figure 6). The bedload particles deflected by the roughness elements further from the wall also impact lower on the wall ( Figure 6a) and impact further downstream on the wall as it takes longer to impact on the wall (Figure 6b).

Calculation of Instantaneous Lateral Erosion Rate
Assuming the channel wall is fully exposed to impacts, the erosion rate E c due to deflections from one grid cell on a roughness element can be expressed as the product of two terms: the volume eroded per particle impact V c and the number of particle impacts per unit time I w (Sklar & Dietrich, 2004): where V c can be calculated as a function of impact velocity v i , and rock parameters, including Young's modulus of elasticity of the bedrock Y, dimensionless bedrock strength coefficient k v , and tensile yield strength σ T : I w can be determined from I c depending on whether the movement of bedload particle deflected by each cell will lead to an impact on the wall or not. If the bedload particle deflected by the roughness element does not impact on the wall, its impact rate on the wall I w is zero. However, if the bedload particle obtains enough momentum to reach the wall, its impact rate on the wall I w is the same with that on the roughness element I c .
E c varies with each grid cell on a roughness element (Figure 7). Only the 1/4 of the semisphere roughness element that faces upstream and toward the near wall contributes to E c due to the concentration of impacts on the upstream facing part of the semisphere ( Figure 5) and the deflection of bedload particles toward the other side of the channel if they impact on the roughness element surface that faces against the wall (Figure 7). E c is highest at the impact position that has a normal vector b n facing 45 degrees relative to the longitudinal centerline of the roughness element in planview and is close to the base of the roughness element, because the rebounded velocity ( Figure 3) and the impact rate ( Figure 5) are highest there. E c decreases with the increasing distance between the roughness element and the wall due to the loss of lateral momentum of bedload particles when traveling toward the wall (Figure 7).
Assuming that bed roughness elements are uniformly distributed in rows composed of equally spaced semispheres (Figure 2), and transported bedload is uniformly distributed across the channel, each row of roughness elements deflects same number of bedload particles and causes same amount of lateral wall erosion. Therefore, only one row of roughness elements is used for calculating the instantaneous local lateral erosion rate E c and the total erosion rate E t due to the existence of one row of roughness elements is simply the sum of all E c , the local erosion rates due to individual bedload particles deflected by each grid cell on the roughness elements: Because the total erosion rate due to multiple rows of roughness elements is the superposition of the lateral erosion rate due to a single row of roughness elements, and the lateral erosion rate in the longitudinal direction repeats for the downstream distance d between two adjacent rows of roughness elements, the integrated lateral erosion rate within d due to multiple rows of roughness element is equal to E t . Therefore, the averaged area of material removed from the channel cross-section per unit time (referred to as bulk erosion rate E b ) within d can be expressed as Bedload particles impact on the wall at many different elevations and downstream locations ( Figure 6). To calculate the average lateral erosion rate E z at a given elevation z, the wall is divided into a uniform grid with a vertical interval Δz from the base of the wall to the maximum erosion height on the wall z lmax . A value of Δz = 1 mm is selected here in accordance with the experimental results of Fuller et al. (2016), and z lmax is obtained from the distribution of the height z l of all impacts.
The impact area within each grid A w is The lateral erosion rate E z for a given elevation range z+Δz can be calculated as a sum of the volume eroded by impacts that fall within that elevation range divided by the impact area A w.

Coevolution of Lateral Erosion Rate, Wall Morphology, and Shear Stress
As the wall is eroded over time, the travel distance, and hence, the potential for loss of momentum of bedload particles after collision with the roughness element, increases, resulting in lower instantaneous lateral erosion rates. Meanwhile, the flow becomes wider and shallower as the wall is eroded. This results in a somewhat lower bed shear stress and hence lower lateral erosion rate. In turn, the lower lateral erosion rate will slow down the wall evolution. Without considering the coevolution between shear stress and lateral erosion rate, the model will exaggerate wall evolution.
To model the effects of wall evolution, we break the simulation into a sequence of time periods, each time period T lasting 10 min. Smaller time periods were tested but did not influence the results. During each period, we assume that the flow depth, and thus shear stress, does not change. We average the erosion rate from impacts that occur during that time period. Then for the next period we update the depth and shear stress and calculate new erosion rates. At beginning of the simulation (T = 1), the initial depth and shear stress are obtained from assuming a rectangular cross-section from equations 1 to 8. As the wall is eroded over time, the channel cross-section and hence the wetted area become irregular. Therefore, Q w (T) is not simply a product of W(T), h(T), and U(T) at the time period T > 1. Instead, Q w (T) needs to calculated from the wetted area A(T) over the irregular cross-section of the flow: where A(T) is a function of flow depth h(T) and needs to be obtained from integrating the flow width over h(T) for a given cross-section shape. We assume that the friction factor f is constant over the run period, because the changes of flow depth are relatively small. Combining equation 36 with equations 1-6, h(T) can be expressed as h(T) and A(T) can be solved from equation 40 by starting with an initial guess of h(T), integrating the flow width over h(T) for the current cross-section shape to get A(T) and iteratively changing the values of h(T) and A(T) until these two solutions converge in equation 37. U(T) will then be back-calculated from equation 39 and used to get the total shear stress τ(T) and shear stress available for sediment transport τ s (T) from total friction factor f = (1 − F r )f b +F r f r and bedrock friction factor (1 − F r )f b using equation 2, respectively: For each period, we calculated the suite of particle deflections and resulting erosion rates, then updated the wall morphology, used it recalculate the water depth, water velocity, and shear stress available for sediment transport in the next time period from equations 36 to 39 and updated the particle impact velocity, particle impact rates, and erosion rates at next period.

Results
We assessed model performance using results from laboratory experiments reported by Fuller et al. (2016). Fuller et al. (2016) constructed three experimental channels (referred to as channels C1, C2, and C3), held the water discharge and sediment supply constant for each channel throughout the experiment, but varied the roughness element size over six classes: no roughness elements (smooth sections); 2.4 mm; 4.3 mm; 7.0 mm; 10.0 mm; and 16.0 mm (roughness sections). Tables 1 and 2 list the initial hydraulic and sediment transport conditions in the Fuller experiments Table 1 Initial Hydraulic and Bedload Transport Conditions Used in the Simulation of the Fuller Experiments One indicates the initial conditions, prior to wall evolution. d Calibrated from τ and τ s by Fuller et al. (2016). e A total of 10.0-mm roughness elements are located both in C2 and C3 by Fuller et al. (2016), the one in C3 is used for calibration of σ T and the one in C2 is used for model performance.

10.1029/2019JF005509
Journal of Geophysical Research: Earth Surface and the values of parameters used in the model calculations. These experiments provide an ideal test case for our model because the flow depth and thus initial shear stress available for sediment transport was measured, and erosion rates and patterns are measured for the various roughness element sizes. However, the rock tensile strength σ T which controls the magnitude of the erosion rate was not measured. For the model calculations we use a value of 5.5 × 10 4 Pa for σ T , which is calibrated from the bulk erosion rate of 10-mm roughness elements (E b = 74 mm 2 /hr) in Channel C3. This value is reasonable for the weak concrete used in the Fuller experiments (Sklar & Dietrich, 2001) and is used for predicting the erosion rate and assessing the model performance for other roughness element sizes.

Model Performance
We assessed three aspects of the model performance when comparing to the Fuller experiments: 1) shape of the eroded profile, 2) peak erosion, and 3) bulk (integrated) cross-section erosion rate. Figure 8 shows the erosion rates measured in the Fuller experiments. An undercut wall morphology occurred, falling within a range of 25-30 mm above the bed for all roughness sections. Lateral erosion was concentrated in the lower half of the undercut (5 to 10 mm above the bed) and decreased progressively up to the maximum height of erosion. The peak erosion was similar for each roughness section, occurring between a height of 0 and 5 mm over 2.15 hr.
The model without coevolving the shear stress, wall morphology, and erosion rate captures the concentration of erosion in the lower half of the wall observed in the Fuller experiments (Figure 8). However, it overpredicts the peak erosion by 3 to 5 times, except the 2.4-mm roughness element where the measured peak erosion is slightly larger (~10%). The lateral erosion concentrates in a smaller zone near the bottom of the wall (below 5 mm above the bed), while the Fuller experiments show a wider concentration zone of erosion spreading from the base of the wall to 10 mm above the bed. The erosion below the radius of bedload particles (4.3 mm) is underpredicted by the model compared with the substantial undercut on the wall observed in the Fuller experiments (Figure 8) due to the assumption of spherical bedload particles, which cannot impact on the wall below their radius.
The Fuller experiments produced a roughly parabolic relation between the roughness element size and integrated cross-section erosion (Figure 9), which increases with roughness size below 4.3 mm, peaks at 4.3 mm, and then gradually decreases with larger roughness element sizes. This relation results from a trade-off between more frequent particle deflections and less shear stress available for sediment transport due to an increase in form drag when the roughness elements get larger (Fuller et al., 2016). Although the model captures this parabolic relation observed in the Fuller experiments, it overpredicts the erosion for all roughness sections by 1.2 to 2 times (Figure 9).
The deviation in the erosion profile and peak erosion rate between the model predictions and the Fuller experiments can occur because changes in wall morphology cause a decline in shear stress applied to the bed. As the wall is eroded over time, the shear stress drops and the travel distance for individual particles increases, resulting in a lower erosion rate over time. We explored the hypothesis that shear stress needs to coevolve with morphology to accurately predict erosion rate by dividing the model run into 10-min periods. Figure 10 shows the decline in mean velocity and shear stress that occurs due to the increase in cross-sectional area as the wall is undercut and the channel is widen. The change in cross-sectional area, velocity, and shear stress is subtle. Shear stress declines most over the first time period (10 min) but barely changes for the rest of the time, because the erosion rate is largest in that first time period, when the bedload particle travel distance is smallest. The overall decline in shear stress is~10%, because the changes in wall morphology are relatively small. Only the bottom of the wall is eroded, and the maximum eroded length is only~10% of the total river width. Our assumption of constant friction factor f may also slow down the change of shear stress.
The model, when coupled with wall evolution, reproduces both the magnitude of erosion and the undercut wall shape that were observed in the Fuller experiments well (Figure 8). The simulated erosion concentrates in the lower half of the erosion zone and tapers off with increasing height on the walls. The predicted peak erosion depth on the wall generally ranges from 8 to 15 mm for all roughness sections, as in the experiments. However, the peak erosion is slightly less than that in the experiments, by~2 mm over the total time period. We suspect this is because we neglected the influence of turbulence on lateral bedload particle deflection into the wall. The Fuller experiments with a planar bed and no deflectors had a wall erosion depth of 2 mm over the 2.15-hr run duration (see Figure 6d by Fuller et al., 2016), while the wall cannot be eroded without deflectors in our model. The model with evolution of the wall and shear stress also successfully reproduces the parabolic relation between the roughness element size and integrated cross-section erosion and the magnitude of erosion over all roughness sections (Figure 9).
The model is suitable for predicting the instantaneous lateral erosion rate on the wall. To successfully predict the change of wall morphology over time, however, the model needs to be coupled with coevolution of shear stress, wall morphology, and lateral erosion rate.

. Evolution of Instantaneous Lateral Erosion Rate and Wall Morphology
The modeled evolution of erosion rate and wall morphology is similar for all channels in the Fuller experiments. Representative profiles, for the 4.3-and 10-mm roughness elements, of lateral erosion rate and wall morphology evolution through time are obtained from our model and shown in Figures 11 and 12, respectively. The instantaneous erosion rate declines over time ( Figure 11). The erosion rate is roughly 10 times lower in the final time period compared to the initial time period. As the wall is eroded over time, the shear stress declines with the mean flow velocity (Figure 9), which leads to a lower erosion rate by decreasing the impact velocity on the wall in later time periods. However, the shear stress at the end of the time period is 90% of the initial shear stress (Figure 9), indicating the influence of the decreasing shear stress on erosion rate is almost negligible over the time period here. The decreasing erosion rate over time is largely due to the

10.1029/2019JF005509
Journal of Geophysical Research: Earth Surface longer travel distance from deflection on the roughness element to impact on the wall as the wall is eroded over time (Figure 12). The erosion rates do not decline to zero over the 2.15-hr model runs, which means the wall can still be eroded if the model continues to run (Figure 12).
At the beginning of the time period, the erosion rate is roughly 10 times smaller in the upper half of the erosion zone, compared to its lower half ( Figure 11). The erosion rate decreases in the lower half of the erosion zone because the undercutting makes the wall further from the roughness elements. In the upper erosion zone, the lower rates of erosion continue because the wall is not as far from the roughness elements at the end of the experimental time. At the end of the time period, the erosion rate in the upper and lower halves of the erosion zone are similar ( Figure 11). The combined effect of this vertical variation through time is a uniform erosion pattern on the wall over the 2.15-hr simulation time ( Figure 12).
The elevation of the peak erosion rate on the wall gets higher from~2.5 to~8 mm above the bed ( Figure 11). Initially, the maximum erosion rate is mostly created by impacts of downward-moving bedload particles, which concentrates in a zone near the base of the wall. As the wall is eroded over time, however, the corner between the bed and the wall is protected as it has been undercut. Instead, more bedload particles will either impact higher on the wall or impact on the bed, obtain upward momentum, and bounce up on the wall. The elevated position of the peak erosion rate on the wall elevates the concentration of erosion zone on the wall (Figure 12).

Coupled Lateral and Vertical Erosion Model
Both field observations (Finnegan & Balco, 2013;Fuller et al., 2009;Hartshorn et al., 2002; and laboratory experiments (Finnegan et al., 2007;Johnson & Whipple, 2010;Shepherd, 1972) have shown that low sediment supply rates promote vertical erosion and high sediment supply rates promote lateral erosion. Vertical erosion is relatively high when bare exposed bedrock is exposed to sediment impact but relatively low when the bed is protected by the alluvial cover. Lateral erosion is thought to be high when the bed is alluviated and able to deflect bedload particles into the wall (Finnegan et al., 2007;Fuller et al., 2016;Gilbert, 1877;Johnson & Whipple, 2010;Shepherd, 1972). However, studies of the competition between lateral and vertical erosion for sediment-flux-driven incision remain qualitative. Our lateral erosion model replicates the essential lateral erosion patterns that were observed in the Fuller experiments by explicitly accounting for bedrock erosion from bedload particle impacts. We couple the lateral erosion model with a vertical erosion model to quantify the changes in vertical and lateral erosion due to impacts from bedload particles for a range of hydraulic and sediment transport conditions. We generalize the lateral erosion model by treating the roughness elements as alluvial cover that has the same grain size as the

10.1029/2019JF005509
Journal of Geophysical Research: Earth Surface bedload particles (D r = D) and use a nondimensional form of the model to show that for a given grain size the full model behavior collapses to a unique functional surface in the parameter space defined by two nondimensional quantities: the relative sediment supply (q s /q t ) and the transport stage (τ * s =τ * c ). We then combine the lateral erosion model with the Sklar and Dietrich (2004) vertical erosion model and quantify the competition between lateral and vertical erosion by looking at the ratio of lateral to vertical erosion rate as a function of relative sediment supply (q s /q t ) and the transport stage (τ * s =τ * c ).

Nondimensional Framework of Coupled Numerical Model
The nondimensional framework for the lateral erosion model is intended to explore the variation of instantaneous lateral erosion rate for the given hydraulic and transport conditions rather than the coevolution of lateral erosion rate, wall morphology, and shear stress over time. We start by determining the size and distribution of roughness elements on the bed. Assuming the alluvial cover provides the only roughness elements capable of deflecting bedload particles and has the same size as the bedload particles (D r = D), the fraction of roughness elements F r increases with sediment supply rate and can be calculated from the relative sediment supply q s /q t using the method proposed by Sklar and Dietrich (2004): where the fraction of roughness elements (alluvial cover). F r is assumed to be a linear function of relative sediment supply. Turowski et al. (2007) developed an exponential formula for F r , and Turowski and Hodge (2017) built a probabilistic framework for the description of the cover effect that contained the linear and exponential models as special cases. For simplicity and consistency with the Sklar and Dietrich (2004) vertical erosion model we elected to use the linear model.
The transport capacity q t can be estimated from the Fernandez Luque and Van Beek (1976) bedload sediment transport relation: where R b = ρ s /ρ w − 1 is nondimensional buoyant density. Assuming the alluvial cover is uniformly distributed on the bed, the distance between two adjacent roughness elements d is expressed as Substituting equation 40 into equation 42, d can be obtained from the given grain size D and relative sediment supply rate q s /q t : We then determine the initial saltation trajectories and deflection trajectories from discrete roughness elements from equations 9 to 15, for a given transport stage τ * s =τ * c and grain size D. These results are then applied in a continuum model by calculating the deflection rates I c on each cell of the roughness surface from equations 16 to 22, the impact rate I w on the wall from I c , the impact velocities v i and positions on the wall from equations 23 to 25 and the resultant total erosion rates E t for all impact locations on the wall from combining equations 29-32 for the given rock parameters (k v , σ T , and Y), relative sediment supply rate q s /q t , transport stage τ * s =τ * c , and grain size D: The downstream velocity after deflection in equation 23 is assumed to be constant here for simplification, without considering the variation of deflection trajectories in the longitudinal direction. To account for the transition from bedload to suspension that is equivalent to a particle taking a hop of infinite length, Sklar and Dietrich (2004) assume that the impact rate on the bed and the impact velocity become negligible as u * approaches w f (see their equation 21 and 22). When l s becomes infinite in our lateral erosion model, the impact velocity on the bed w s (equation 12) before deflection and, hence, the impact velocity on the wall v i (equation 30) monotonically increases with higher transport stage. This is problematic because the lateral erosion rate should decline as the transport stage approaches the suspension threshold. To keep the lateral erosion model consistent with the Sklar and Dietrich (2004) vertical erosion model, v i is set to be negligible by multiplying it with (1 − (u * /w f ) 2 ) 0.5 in equation 44 as u * approaches w f and rearranging equation 44: To evaluate the average lateral erosion rate E l on the wall, E t is averaged over the maximum impact elevation on the wall z lmax which is obtained from the distribution of z l of all deflection trajectories on the wall: The variable u * /w f is a function of transport stage τ * s =τ * c for a given grain size D, so E l is a function of four variables, including rock parameters (σ T and Y), relative sediment supply rate q s /q t , transport stage τ * g =τ * c , and grain size D. The influence of rock parameters (σ T and Y) in equation 46 can be erased when E l is nondimensionalized as (Sklar & Dietrich, 2004) Therefore, E * l can be considered as a function of just two nondimensional quantities, the relative sediment supply q s /q t and the transport stage τ * s =τ * c for a constant grain size D. Meanwhile, an analytical solution for the nondimensional vertical erosion rate E * v has been proposed to be a function of q s /q t and τ * s =τ * c by Sklar and Dietrich (2004): Vertical and lateral erosion can be coupled from the ratio e: because both erosion rates can be related to two variables q S /q t and τ * s =τ * c for a given D.

Competition Between Vertical and Lateral Erosion
In order to explore the competition between vertical and lateral erosion with varied q S /q t and τ * s =τ * c , we assume that channel erosion is disconnected from the hillslopes. The most direct analog for the coupled model here is a bedrock canyon or gorge that is deeply incised into a river valley and largely disconnected from the hillslopes. In order to implement lateral and vertical erosion in a coupled format, we must specify various parameters, including the grain size of transported material, transport thresholds and various sediment, and rock and water properties. For convenience, we use values reported by Sklar and Dietrich (2004) for the South Fork Eel River in Northern California (Table 3).
The first step in exploring the competition between vertical and lateral erosion involved calculating how E * v varies with q s /q t and τ * s =τ * c for the grain size D = 0.06 m at the reference site, using the Sklar and Dietrich (2004) model. E * v has an analytical solution for q s /q t and τ * s =τ * c , so we can simply determine E * v for each combination of q s /q t and τ * g =τ * c from equation 45. Figure 13 shows that E * v collapses to a unique

10.1029/2019JF005509
Journal of Geophysical Research: Earth Surface functional surface in the parameter space created by q s /q t and τ * s =τ * c . As in Sklar and Dietrich (2004), E * v goes to zero at the threshold of motion and suspension along the τ * s =τ * c axis and the threshold of full cover and no cover along the q s /q t axis. The decline in erosion rate at the threshold for suspension is adopted for simplicity here, but we recognize that this is not strictly correct and that there is some reduced bedrock erosion beyond the suspension threshold (Lamb et al., 2008;Scheingross et al., 2014). E * v peaks at the intermediate transport stages (Figure 14a) where the growth in the impact energy is balanced by a decline in the impact frequency as the saltation hop length increases with shear stress, and at moderate relative sediment supply (Figure 14b), where the growth in impact rate is balanced by the reduction in the extent of bedrock exposure with increasing sediment supply.
The second step in examining the competition between vertical and lateral erosion was to explore how E * l varies with q s /q t and τ * s =τ * c for the grain size D = 0.06 m at the reference site. We varied τ * s =τ * c from 1 to 22, and for each value of τ * s =τ * c calculated the initial saltation trajectories (equations 9-12) before deflection by roughness elements and the transport capacity q t (equation 38). We also varied q s /q t from 0 to 1 and for each value of q s /q t calculated the distance between two adjacent roughness elements d (equation 40). For each combination of q s /q t and τ * s =τ * c , we calculated the sediment supply rate q s (equation 37) and used the deflection model to get all the possible individual deflection trajectories from discrete parts of the roughness elements (equations 13-15). We then applied these results in the continuum model by calculating the deflection Using this nondimensional framework, the lateral erosion model also collapses to the unique functional surface in the parameter space defined by q s /q t and τ * s =τ * c (Figure 15). Figure 15 reveals that E * l goes to zero at the threshold of motion and suspension along the τ * s =τ * c axis, and the threshold of no cover along the q s /q t axis, but is relatively high at the threshold of full cover. As with E * v , E * l peaks at an intermediate transport stages, however, E * l peaks at high relative sediment supply rate (~0.7; Figure 15). Figures 14c-d illustrates the pattern of E * l with increasing shear stress and relative sediment supply rate more clearly. E * l shows a parabolic variation with transport stage, where E * l is zero at the threshold of motion due to a lack of particle movement along the transport stage axis (Figure 14c). As the transport stage exceeds the threshold for motion, E * l increases gradually with transport stage due to the growth in impact velocity. However, the impact frequency of bedload particles on the roughness element decreases with transport stage, because the saltation trajectories tend to grow more elongated with increasing shear stress. The growth in the particle impact energy and the reduction in the impact frequency with increasing shear stress results in a peak E * l at intermediate transport stages. E * l goes to zero at the threshold of suspension, because no impacts between roughness elements and bedload particles occur in our model. This is an artifact of the saltation model used. Some limited lateral erosion is possible from deflected particles above the suspension threshold, but E * l would be low. Along the relative sediment supply axis, a parabolic variation of E * l also exists. E * l is zero when the bed is free of cover and remains negligible when the relative sediment supply is <0.15 Figure 14. Nondimensional vertical erosion rate as a function of (a) transport stage τ * s =τ * c and (b) relative sediment supply q s /q t ; nondimensional lateral erosion rate as a function of (c) transport stage τ * s =τ * c and (d) relative sediment supply q s /q t and the ratio of lateral to vertical erosion rate as a function of (e) transport stage τ * s =τ * c and (f) relative sediment supply q s /q t . (Figure 14d). This occurs because when the relative sediment supply is low, the fraction of bed coverage is low, and there are relatively few deflectors on the bed. E * l gradually grows with the relative sediment supply rate above 0.15 due to the increase of the number of saltating bedload particles and the extent of roughness. However, E * l peaks at the relative sediment supply of~0.7 (Figure 14d) due to a competition between the impact area A c and wall-normal velocity v o and the number of deflections on each cell of the roughness surface. When q s /q t increases, the distance between two adjacent roughness elements starts to decline, which will reduce the deflections near the bottom of the roughness surface and force bedload particles to impact near the top of the roughness surface. The concentration of impacts near the top of the roughness surface will lead to lower impact area on each cell as the cell starts to get close to the flux surface ( Figure 5) and lower wall-normal velocity v o after deflection by the cell as the vertical velocity w s before deflection declines and the normal vector increasingly points upward. However, the number of deflections on each cell increases with higher sediment supply rates as q s /q t increases. The decrease in impact area A c and wall-normal velocity v o and the increase of the number of deflections on each cell of the roughness surface with higher q s /q t will lead to a peak in E * l when they are balanced. E * l starts to decline for q s /q t above~0.7 and is~75% of the peak lateral erosion rate at the threshold of full cover.
The contour lines of nondimensional lateral erosion rate are not smooth. This is not improved by using smaller time steps and space grids. The roughness element surface is discretized into nearly uniform triangular grid cells to model the collision with a finite number of bedload particles, which leads to variations in the modeled erosion rate. Some variation is also caused by our numerical approach. We track the movement Figure 15. Nondimensional lateral erosion rate (E * l ) as a function of transport stage and relative sediment supply.
of each particle to obtain the impact velocity and the impact position on the wall under every combination of relative sediment supply rate and transport stage instead of deriving explicit empirical correlations, resulting in a lateral erosion model that varies irregularly with control variables.
The competition between vertical and lateral erosion was calculated from the ratio of E * l to E * v for each combination of q s /q t and τ * s =τ * c . The ratio e collapses to a unique functional surface in the parameter space created by q s /q t and τ * s =τ * c (Figure 16). e goes to zero with no bed cover, at the thresholds of motion and suspension, and is infinite when the bed has full cover. Figures 14e-f illustrates the patterns in e with changes of τ * s =τ * c and q s /q t . Along the τ * s =τ * c axis, e is parabolic, with a peak at an intermediate transport stage. This occurs because E * l increases more rapidly than E * v at lower transport stages (<10) and decreases more rapidly than E * v at high transport stages (Figures 14a-b), for a constant q s /q t . In contrast, e shows a monotonic increase with increasing q s /q t (Figure 14f); e goes to zero when q s /q t = 0 and gradually increases with relative q s /q t (>0.15), because E * l grows faster than E * v when the relative sediment supply rate is below 0.5, and E * l continues to increase but E * v starts to decrease when the relative sediment supply is between 0.5 and 0.7 (Figures 14b  and 14d). The ratio e continues to increase at high relative sediment supply (>0.7), because E * l decreases more slowly than E * v . When the bed is fully covered, the ratio goes to infinity as the lateral erosion rate is relatively high, but the vertical erosion rate goes to zero.
The coupled model shows that the lateral erosion rate is lower than the vertical erosion rate under nearly 75% of the transport and supply conditions (Figure 16). Lateral erosion is negligible at low sediment Figure 16. The ratio of lateral to vertical erosion rate e¼E * l =E * v as a function of transport stage and relative sediment supply.

10.1029/2019JF005509
Journal of Geophysical Research: Earth Surface supply rates when the bed coverage is less than 20% and gradually increases with the extent of alluvial cover but only dominates at high sediment supply rates when the bed is largely covered by roughness elements.
The ratio e is ultimately controlled by the change in E * v and E * l and where it is high does not necessarily correspond to where either E * v or E * l are largest. Nevertheless, lateral erosion only dominates over vertical erosion under a limited range of conditions.

Discussion
The lateral erosion model confirms that bedload particle impact is a viable mechanism for lateral erosion in bedrock rivers by reproducing key patterns in lateral erosion from the Fuller experiments, including the undercut wall shape, the peak erosion, and the total erosion rate. Saltating bedload particles obtain lateral momentum to erode the wall by colliding with the roughness elements on the bed. The bedload particle impacts concentrate in a zone near the bottom of the walls, thereby creating an undercut wall shape.

Limiting Conditions on Lateral Erosion
While our model can reproduce key features of lateral rock erosion in channels, it is useful to consider some limiting conditions on the process of lateral erosion by abrasion. During periods when wall erosion can be effective, there are limits to how far lateral erosion by abrasion may occur before one of the following happens: 1) changes in channel geometry cause the stress to fall below the threshold of motion to maintain bedload; 2) the undercut becomes so deep that deflected particles can no longer reach the wall; or 3) the undercut is so deep that the rock mass above it fails into the channel (as in Figure 1a). It is unlikely that the first limitation will actually occur in a river undergoing undercutting. As the wall is undercut over time, mean velocity and shear stress drop due to the increase in cross-sectional area. The lateral erosion rate will go to zero when the shear stress is below the threshold for particle motion. However, this is unlikely to happen because the stress and wall morphology coevolve. At low stresses, where changes in the wall would affect the shear stress, the erosion rate would be low, so the wall evolution would be very slow. It would therefore take an excessively long time for the shear stress to drop below the threshold of motion.
The lateral momentum for bedload particles to reach the wall drops due to the increase in travel distance as the wall gets undercut over time. To explore how this affects further undercutting, we ran the lateral erosion model over 15 hr using the 10-mm roughness element section. The lateral erosion stops after 12 hr, although the transport stage (~2.5) at the end of the time period is still enough to transport bedload particles (Figure 17). This occurs because the wall is eroded over 18 mm at the

10.1029/2019JF005509
Journal of Geophysical Research: Earth Surface end of the time period, which is too far for bedload particles to impact on the wall at transport stage of 2.5. As such, the increase in travel distance provides a greater limiting condition on lateral erosion than the drop of shear stress. Using a constant resistance coefficient over the run may overpredict the shear stress because the hydraulic roughness may increase as the wall is undercut. Also, the calculation shown in Figure 17 does not include any roughness elements within the undercut. If the newly exposed bed by channel undercutting becomes alluviated, those deflectors would allow lateral erosion to continue.
Continued undercutting of the lower wall creates an imbalance on the wall and may cause the upper part to collapse and to widen the whole channel. Such a mechanism of channel widening has been documented in both experiments (Carter & Anderson, 2006) and field data (Cook et al., 2014). However, the question of how far the wall needs to be undercut before it fails remains unanswered. Bedrock walls with lesser rock mass strength can fail more easily as the lower part of the wall is undercut. The degree of fracturing and jointing on the bedrock walls influences the rate of rock sliding and toppling and hence channel width. Bedrock bedding may play a dominant role in controlling the wall collapse. Undercut bedrock walls with vertical bedding can cause a channel to widen more effectively than with horizontal bedding, which may remain intact for deeper undercuts.

Undercut Wall Shape Dynamics
One of the key findings of our model is that in bedrock channels with a planar bed, the competition between vertical and lateral erosion is controlled by the extent of alluvial cover under different sediment supply conditions, which agrees with the experimental and field measurements (Finnegan et al., 2007;Finnegan & Balco, 2013;Fuller et al., 2009Fuller et al., , 2016Johnson & Whipple, 2010). Our model provides the opportunity to quantify how different wall shapes are formed under different sediment supply conditions. In a low sediment supply environment, the channel bed is more exposed and vertical erosion will dominate, with lateral erosion relatively negligible, resulting in a near straight wall shape. At an intermediate to high sediment supply where the bed is 50%-90% covered, both the bed and walls can be cut by bedload particle impacts. The continuing lowering of the channel bed will shift down the lateral erosion zone by preventing the bedload particles impacting on a fixed elevation on the walls. This will create an undercut wall shape that keeps the same width but spreads more deeply over time. However, when the bed is near fully covered (>90%), the bed is relatively static due to the protection of alluvium, leading to an undercut wall shape that gets wider over time. As such, the wall shape would change from near straight to deeply undercut as the sediment supply increases.
The undercut wall shape may be modified by roughness elements made of the bedrock surface. The bedrock rivers have a wide range of sculpted bed morphologies (Montgomery & Buffington, 1997;Richardson & Carling, 2005;Wohl, 1993;Wohl & Merritt, 2001), such as potholes, flutes, furrows, runnels, etc. In a bedrock channel with bedrock obstacles near the walls, bedload particles can be deflected toward the walls by bedrock obstacles even when no alluvial cover exists. Beer et al. (2017) mapped the lateral erosion patterns in a bedrock gorge in the Swiss Alps under three bedrock obstacle conditions: 1) no bedrock obstacle; 2) low bedrock obstacle; and 3) high bedrock obstacle. Bedrock walls without bedrock obstacles were only slightly eroded, while in sections with low and high bedrock obstacles, the walls were undercut deeper and higher above the bed. The occurrence of bedrock obstacles to deflect bedload particles to higher elevations than the alluvium may have the effect of elevating the undercut zone. The size of bedrock roughness obstacles can influence the erosion rate from two opposite effects. Small bedrock obstacles do not have large surface area to deflect bedload particles but tend to have high impact velocity due to low form drag. Larger bedrock obstacles have more surface area for deflections, but the impact velocity will be reduced because of higher form drag. Intermediate bedrock obstacles that balance the trade-off between surface area and impact velocity may result in highest lateral erosion rates.
One implication of our model would suggest is that the extreme magnitude of supply events may not be the most effective in eroding bedrock walls. Instead, the more frequent high moderate sediment supply events may be responsible for the majority of bedrock wall erosion. This occurs because the reduced impact area on each roughness element and the reduced impact velocity with increasing alluvial cover (>75% coverage; Figure 15).

Journal of Geophysical Research: Earth Surface
It is possible to infer the relative width to depth ratio and degree of incision of a channel cross-section from Figure 16. A bedrock channel with a high sediment supply rate, which can be found near the upper part of Figure 16, is mostly covered by alluvium. This channel would be dominated by lateral erosion with negligible vertical erosion, allowing for a wide bedrock channel, relative to its depth. In contrast, a channel that receives relatively little sediment supply should plot near the lower part of Figure 16, will preferentially incise the bed, and have a lower relative width to depth ratio. Of course, the sediment supply and transport stage conditions of bedrock rivers change over time with hydrographs and sedigraphs in a basin. The ultimate shape of a channel is determined by how long it spends in particular positions on Figure 16. A channel that spends the vast majority of its time in the lower corner of Figure 16 is likely to be narrow and deeply incised. A channel that is in the upper corner of Figure 16 most of the time will be relatively wider. Tracking a channel through time on Figure 16 requires a full morphodynamic implementation of the model presented herein, which requires imposed hydrographs and sedigraphs.

Model Limitations and Further Prospects
There are a number of simplifications in our model, which were necessary to produce a result, that may affect the outcomes. Our model is a sediment flux-driven incision model that is most applicable to rivers flowing through massive crystalline rock that does not take into account any other lateral erosion mechanisms. Plucking may dominate bedrock erosion where rock is weak or well jointed (Chatanantavet & Parker, 2009 or where shear stresses are large enough to entrain or transport large blocks (Beer et al., 2017;Lamb et al., 2015;Montgomery, 2004 ;Stock et al., 2005). Plucking is particularly effective in erosion of columnar basalts (Dubinski & Wohl, 2013;Larsen & Lamb, 2016c) and may be possible by fluid shear stresses plucking pieces of weak rocks from channel walls (Montgomery, 2004;Stock et al., 2005).
Our model uses a uniform grain size with spherical shape for sediment particles to represent the wide distribution of grain sizes supplied to bedrock rivers. Grain size controls the threshold for motion and hence the transport stage and, hence, impact velocity and impact rates. Grain size of the alluvial cover determines the elevation of collision, thereby influencing the transfer of momentum during collision and the impact height on the wall. High points of the alluvial cover that protrude above the saltation layer do not deflect bedload particles into the wall. Therefore, the distribution of grain sizes supplied by the upstream catchment (Sklar et al., 2017) may influence the lateral erosion rate by changing the fraction of total load that is transported as bedload and the momentum transfer of bedload particles during collision with the alluvial cover. The shape of sediment particles determines the distribution of impact angles during collision between roughness elements and bedload particles, thereby influencing the direction of movement after collision. Given that our assumption of a uniform grain size with spherical shape has well reproduced the erosion patterns observed in the Fuller experiments, which used nonspherical deflectors, the influence of the nonspherical shapes of natural particles on lateral erosion rate may be negligible. Our lateral erosion model uses numerical formulations to track the movement of individual bedload particles. The potential for bedrock erosion by bedload impacts at transport stages above the suspension threshold is ignored. It is possible that particle impacts might be viscously damped for fine grains that are transported as suspended load. Yet bedload transport remains a significant but decreasingly important component of the total load as transport stages increase above the suspension threshold (Lamb et al., 2008;Scheingross et al., 2014). Suspended load has been proposed to be responsible for lateral erosion through turbulent fluctuations that laterally sweep particles to impact on the wall (Whipple et al., 2000). It is not possible for us to track particle movements above the suspension threshold, so we force the lateral erosion rate to zero at the suspension threshold, which is consistent with the Sklar and Dietrich (2004) vertical erosion model. Further work is necessary to develop a version of our model to handle lateral erosion by sediment suspension.
The simplified treatment of flow dynamics in the model may influence the result as well. Movement of sediment after collision is modeled by assuming that the influence of turbulence on trajectories is negligible. However, local turbulent fluctuations can be intense above a bed with significant roughness (Richardson & Carling, 2005). We assume that flow advection is negligible near the bed so that particles impact on roughness elements and subsequently on the wall without being swept away with the flow. The advective component of the impact velocity can be significant over roughness elements (Johnson & Whipple, 2007;Tinkler, 1997), where flow goes around large roughness elements and advects the sediment toward the wall, potentially increasing the impact velocity and rates on the wall. Some caution should be exercised in applying the model where the cross-section is irregular or where the flow field is nonuniform.
An important simplification in the lateral erosion model is our assumption that bedload particles are uniformly transported in a rectangular channel with a planar bed and straight walls. At low relative sediment supply rates, the bedload layer self organizes into a concentrated filament along the channel centerline. In a rectangular bedrock channel, the shear stress is higher in the channel center than near the wall due to the wall drag (Parker, 1978). This flow structure results in faster bedload particle velocity in the channel center than near the walls. Nelson and Seminara (2011) modeled the Finnegan et al. (2007) experiments and showed that the formation of a bedload filament along local channel depressions is capable of forming an incised channel with strath terraces because bedload is gravitationally drawn to the lowest part of the channel cross-section, forming preferential pathways for bedload in bedrock channels (Chatanantavet & Parker, 2008;Finnegan et al., 2007;Inoue et al., 2014;Inoue et al., 2016;Nelson & Seminara, 2011Turowski & Hodge, 2017). The higher speed and greater concentration of bedload particles in the channel center will increase the impact energy and frequency and accelerate the vertical erosion rate in the channel center but slow down the lateral erosion rate due to the increasing travel distance for the particles to impact on the wall. This can have a self-reinforcing effect where the bedload pathways are topographically steered by an incised groove so that bedload concentrates vertical erosion in the channel center until the incised groove is alluviated (Cao, 2018). Alternatively, the areal concentration of bedload can enhance the lateral erosion on the groove walls because sediment particles can be deflected by the alluvium to impact on the sidewalls and to widen the incised groove. The organization of bedload layers is not presently considered in our model because there is no analytical solution to predict this phenomenon, but this needs to be addressed to understand complex channel cross-sections.
Our model also does not consider the self-organization of bedload layers that may lead to the formation of persistent alluvial patches at high relative sediment supply rates. Relative sediment supply in the Fuller experiments was kept below the patch formation threshold. Field observations and flume experiments have shown that partial alluvial cover tends to self-form in patches in bedrock rivers and may form alternating gravel bars (e.g., Chatanantavet & Parker, 2008;Lisle, 1986). Turowski (2020) developed a lateral erosion model that more explicitly treats the effects of self-organized alluvium. The model does not track particle deflections and evolution of the wall, as in our model. Instead, a deflection length is defined as saltation hop length times its defection angle and if the deflection length is less than the distance between the bar and the opposite wall, lateral erosion is thought to occur, increasing the channel width to a steady state reach width. Our model would predict the opposite erosion pattern predicted by the Turowski (2020) model. Our model would underpredict erosion on the wall adjacent to a bar because the bar will have densely packed deflectors, and it would overpredict erosion on the wall opposite to the bar because the bare bedrock will have few immobile deflectors. Inclusion of flow patterns associated with alternate bars may produce the erosion patterns predicted by the Turowski (2020) model. Flow in bends directs coarser particles toward the outer bank (Dietrich & Smith, 1984;Dietrich & Whiting, 1989), which has been shown experimentally to enhance lateral erosion downstream of the bend apex (e.g., Mishra et al., 2018). Future development of our model will need to more faithfully treat the self-organization of alluvial cover and flow in bends to predict the effects of channel curvature and self-organized alluvium.
Despite the simplifications in our model, it agrees well with experiments that have relatively simple geometries (Fuller et al., 2016), which suggests that the model captures the fundamental mechanism correctly. Furthermore, the model generates undercut wall shapes that are qualitatively similar to field observations (Beer et al., 2017). Application of the model to a natural channel needs to consider timescales of effectiveness for both the vertical and lateral erosion processes, which are ultimately controlled by discharge and sediment supply variations (Finnegan et al., 2005;Finnegan & Balco, 2013;Inoue et al., 2014Inoue et al., , 2016Lague, 2010;Lague et al., 2005). Wall erosion is the integrated result of intermittent periods of variable discharge and sediment supply. Ultimately, application of our model to natural channels with variability in discharge and sediment supply requires a morphodynamic model by 1) developing a lateral erosion model by suspended load; 2) finding an analytical solution of the model to make the calculations tractable through geologic time; and 3) developing a way to parameterize our flow resistance submodel.

Conclusion
We have developed a mechanistic model for lateral erosion of bedrock channel banks by bedload particle impacts using well-established empirical relations for initial velocities of bedload particles, a simplified reflection methodology for collision with roughness elements, and a numerical model for tracking the motion of bedload particles from collision to impacts on the wall. Simulations of the Fuller experiments show that the model successfully predicts the essential undercut wall shape, the dynamics of peak erosion rate and total cross-sectional erosion rate with roughness element size, which not only validates the formulation of our lateral erosion model but also supports the bedload particle impacts as an effective mechanism for lateral incision in bedrock rivers. The predicted lateral erosion rate can be further expressed in nondimensional form as a function of transport stage and relative sediment supply for the given grain size by assuming that the alluvial cover due to deposition of sediment particles is effective at deflecting downstream transport particles. The nondimensional lateral erosion model defines a unique functional surface bounded by four thresholds, including the threshold of motion, the threshold of suspension, the threshold of no cover, and the threshold of full cover. The lateral erosion is relatively high at the threshold of full cover but turns to be zero at all other three thresholds. The model also predicts a peak lateral erosion rate when the bed is near 70% covered, due to a trade-off of deflection rates and deflection angles as the sediment supply increases. A coupled model that combines vertical erosion with lateral erosion due to bedload particle impacts is further developed. The coupled model predicts that vertical erosion dominates under~75% of transport and supply conditions on the unique functional surface. The lateral erosion only outpaces the vertical erosion when the bed is near fully covered.

Notation
A c projected area of gid cell (m 2 ) A w impact area on the wall (m 2 ) A(T) wetted area at time period T (m 2 ) C 1 gravitational acceleration coefficient (ms −2 ) C 2 drag deceleration coefficient (ms −2 ) C d drag coefficient (dimensionless) C r restitution coefficient (dimensionless) d distance between two adjacent roughness elements (m) D grain size of bedload particles (m) D r grain size of roughness elements (m) e the ratio of lateral to vertical erosion rate (dimensionless) E c instantaneous total erosion rate from one cell (m 3 s −1 ) E b bulk erosion rate (m 2 s −1 ) E l lateral erosion rate (ms −1 ) E * l nondimensional lateral erosion rate (ms −1 ) E t total erosion rate (m 3 s −1 ) E * v nondimensional vertical erosion rate (ms −1 ) E z erosion rate at elevation z (ms −1 ) f total friction factor (dimensionless) f a friction factor for alluvium (dimensionless) f b friction factor for bedrock surface (dimensionless) f r friction factor for roughness element (dimensionless) F a fraction of alluvium (dimensionless) F r fraction of roughness elements (dimensionless) g gravity acceleration (ms −2 ) h water depth h c impact elevation on the roughness element (m) h s saltation hop height (m) I c impact rate on grid cell (s −1 ) I p impact rate on the upward trajectory plane (m −2 s −1 ) I w impact rate on the wall (s −1 ) i s incoming saltation velocity vector (ms −1 ) IL impact position vector (ms −1 ) IV impact velocity vector (ms −1 ) k s hydraulic roughness length scale (m) k v Rock resistance coefficient (dimensionless) l s saltation hop length (m) l sd saltation hop length of the downward trajectory (m) l su saltation hop length of the upward trajectory (m) L length of upward trajectory plane (m) M mass of bedload particles (kg) b n surface normal vector (dimensionless) N number of grids (dimensionless) o s outgoing saltation velocity vector (ms −1 ) p projection of the incoming velocity vector onto the surface normal vector (ms −1 ) P(T) wetted perimeter at time period T (m) q s sediment supply per unit width (kgm −1 s −1 ) q t transport capacity per unit width (kgm −1 s −1 ) Q w water discharge (m 3 s −1 ) r semi-circle radius cut along the roughness element (m) R b nondimensional buoyant density (dimensionless) S channel slope (dimensionless) S t Stokes number (dimensionless) u instantaneous downstream velocity (ms −1 ) u * shear velocity (ms −1 ) u i downstream impact velocity (ms −1 ) u o outgoing downstream velocity (ms −1 ) u s saltation downstream velocity (ms −1 ) U z flow velocity at depth z (ms −1 ) U(T) flow velocity at time period T (ms −1 ) v instantaneous lateral velocity (ms −1 ) v min minimum wall-normal velocity limit (ms −1 ) v i lateral impact velocity (ms −1 ) v o outgoing lateral velocity (ms −1 ) v s saltation lateral velocity (ms −1 ) V l volume eroded per particle impact (m 3 ) w instantaneous vertical velocity (ms −1 ) w 0 outgoing vertical velocity (ms −1 ) w f fall velocity (ms −1 ) w i vertical impact velocity (ms −1 ) w s saltation vertical velocity (ms −1 ) W channel width (m) x l downstream impact position (ms −1 ) y l lateral impact position (m) Y Young's modulus of elasticity (kgm −1 s −2 ) z l vertical impact position (m) z lmax maximum erosion height on the wall (m) ρ s sediment density (kgm −3 ) ρ w water density (kgm −3 ) σ T Rock tensile strength (Pa) τ total shear stress (Pa) τ s shear stress due to skin friction (Pa) τ * s nondimensional shear stress due to skin friction (dimensionless)