Numerical Simulations of Meanders Migrating Laterally as They Incise Into Bedrock

The problem of meandering in mixed bedrock‐alluvial rivers is more challenging than that of purely alluvial streams, in that alluvial, bed incisional and bank incisional morphodynamics must be accounted for. Here we present a numerical formulation that addresses heretofore unanswered questions. Bed incision is based on abrasion due to saltating grains. The model satisfies mass conservation of alluvium over a partially covered bedrock surface. Bank incision is treated in terms of a measure of the incipient collision of bedload particles with the bank. It is assumed that land accretes along the inside of point bars when the water depth falls below a specified shallow value. All but one of the runs are performed with repetitive two‐step hydrographs. Runs starting from a low‐amplitude sine‐generated curve indicate that sinuosities at least as high as 2.5 can be achieved. The rate of increase of sinuosity declines in time, but does not vanish. For the same hydrograph, increasing the initial thickness of alluvium on the bed causes the rate of vertical bedrock incision to decline, and bend sinuosity to increase at a faster rate. At a sufficiently high thickness, the channel migrates laterally without bed lowering. Bend shape can change dramatically with increasing alluvial thickness, with high thickness favoring more regular bend trains. For the same initial alluvial thickness, increasing the peak flow of the hydrograph causes the vertical incision rate and the rate of sinuosity growth to increase. The model thus captures a wide range of behavior associated with bedrock meandering.

dissipation (Langston & Tucker, 2018). So far, two models have been proposed for lateral bedrock erosion caused by the collision of bedload particles. Fuller et al. (2016) modeled the relationship between the lateral erosion rate and the roughness of the bed surface, because bedload particles moving in the streamwise direction are irregularly deflected into the wall by bed and wall roughness elements. Their experimental results indicated that the lateral erosion rate increases with increasing bed roughness. Recently, Li et al. (2020) proposed a numerical model for tracking the motion of bedload particles that are deflected by roughness elements and collide with walls, and succeeded in reproducing the experiment by Fuller et al. (2016). Inoue, Parker, et al. (2017) modeled the lateral erosion rate using the lateral bedload transport generated by flow meander and bar development. Field observations of the Daan gorge (Taiwan) by Cook et al. (2014) showed that the lateral erosion rate in a straight section of a channel is smaller than that found in an adjacent channel bend. Mishra et al. (2018) experimentally showed a tendency similar to that of Cook et al.'s observations and revealed that the development of a point bar has a great influence on the lateral erosion in the channel bend. These observations generally support the concept of the model proposed by Inoue, Parker, et al. (2017).
The scalar value of the lateral bedload transport rate is modeled as zero at the bank in an Eulerian dispersion model (Hasegawa, 1989;Mosselman & Crosato, 1991) because the number of colliding particles and the number of rebounding particles are equal. Therefore, Inoue, Parker, et al. (2017) proposed a model in which the lateral erosion rate y E of a bedrock sidewall is dependent on the near-wall gradient of the lateral bedload transport rate; where by q is the bedload transport rate in the lateral direction,  bank is a nondimensional abrasion coefficient for bedrock sidewalls, and b L is an estimate of the thickness of the boundary layer over which the transverse bedload rate drops to zero at the bank ( Figure 2). Secondary flow and lateral bed slope determine the lateral bedload transport rate in channel bends with point bar formation (Hasegawa, 1989  to, 1991). Secondary flow near the channel bend forces the near-surface fluid toward the outside of the bend, and the near-bed fluid toward the inside point bar, thus, pushing the bedload toward the inner bank. On the other hand, the lateral bed slope pushes sediment toward the outer bank. The lateral sediment transport rate is zero in an equilibrium state where the influence of the secondary flow and the influence of the lateral bed slope are fully balanced. This suggests that the lateral erosion of bedrock sidewalls may eventually become zero even in a meandering channel. There may be other factors that direct the particles toward the side walls, for example, fluctuations in flow discharge, lag due to fluid inertia, irregular motion of bedload particles, etc.
Observing the evolution of incising bedrock meanders in the field is daunting, as they evolve over a time scale of hundreds to millions of years (Hancock & Anderson, 2002). As a result, numerical methods provide an efficient tool for investigating the factors responsible for determining the morphology of meandering in bedrock (Limaye & Lamb, 2014).
Although several numerical models have been proposed to study the evolution of bedrock meanders (Howard & Knutson, 1984;Limaye & Lamb, 2014, these models do not explicitly treat river channel sediment transport and impact wear of bedrock sidewalls. Inoue, Parker, et al. (2017) presented a model describing the morphodynamics of a mixed bedrock-alluvial river bend, including sediment transport. Their simulations demonstrated the existence of an approximate state of dynamic equilibrium for the cross-sectional shape of a channel with streamwise constant curvature (uniform bend). In their simulations, however, the migration rate of the outer bank is a specified parameter. The lateral migration rate should, more properly depend on the sediment transport rate, channel curvature, bedrock bank strength, etc.
In this study, we implement the bedload impact abrasion model for a bedrock bank embodied in Equation 1 into the two-dimensional morphodynamics model proposed by Inoue, Parker, et al. (2017). Then, we investigate the lateral sediment transport rate in the dynamic equilibrium state of the riverbed, as well as the effects of sediment supply and flow discharge on the shape of meanders incised into bedrock. Finally, we discuss several aspects of meander morphology (e.g., terrace formation, channel sinuosity, and skewness) in mixed bedrock-alluvial channels.
Although we do not try to reproduce any particular river reach using our model, Figure 1 illustrates natural morphologies of meandering bedrock rivers that have inspired our model. Figure 1b shows the meandering bedrock-alluvial Shimanto River. Lateral migration is documented by the many meander cutoffs along the river; these have become perched due to uplift, and are used for farmland or towns. The hillslope adjacent to an inner bank gradually falls into the river, but in the interim can be used for farmland or towns. Outer banks are not overhanging, but are steeper than inner banks. Figure 1c shows a tributary of Sugar Creek, Indiana. Evidently, the sidewalls are tougher here, forming overhangs as a bend migrates into them.

Numerical Model
The equations governing the flow field are the two-dimensional depth-averaged equations for mass and momentum conservation (see Appendix A). Incisional/alluvial morphodynamics is calculated by employing the framework proposed by Inoue et al. (2014) (see Sections 2.1 and 2.2). The differences between our model and the model proposed by Inoue, Parker, et al. (2017) pertain to the lateral erosion model presented in Section 2.3 and the land accretion model presented in Section 2.4.
Our study uses equations in a moving boundary-fitted coordinate (MBFC) system (Asahi et al., 2013), but for simplicity, we express the equations here in an orthogonal coordinate system (x-y coordinate system). The effects of secondary flow are not included in the depth-averaged flow simulation, but are included as vector changes in the bed load transport rate (see Equations 11 and 12). As the equations for the bedload vector have been proposed in the  s n coordinate system, where s is the streamwise direction and n is the normal direction, we show in Equations 4 and 5 how to convert them to the  x y coordinate system of the computational grid.

Alluvial Bed Deformation
To calculate the thickness of the alluvial layer above bedrock, we use a form of sediment conservation that specifically includes sediment moving in the bedload layer, in addition to sediment deposited on the bed (Inoue et al., 2014, Inoue, Parker, et al., 2017: V is the volume moving as bedload per unit area,  is the porosity of the alluvium, c p is the areal fraction of alluvial cover,  a is the thickness of the alluvial layer, and bx q and by q are the bedload transport rates per unit width in the x and y direction, respectively. In alluvial rivers, b V is equal to the saturation volume of bedload per unit area, bs V (i.e., the volume at bedload transport capacity). In a mixed alluvial-bedrock channel, however, b V is not necessarily equal to bs V . Hence, bx q and by q are parameterized as follows:  where bcx q and bcy q are the capacity bedload transport rates in the x and y directions, respectively.
These capacity components are calculated using the relations: where bcs q and bcn q represent the bedload transport capacity per unit width in the streamwise and transverse directions respectively, and x U and y U are depth-averaged velocities in the x and y directions obtained from the calculation of the two-dimensional flow field.
The streamwise bedload flux bcs q can be estimated using the empirical relation of Meyer-Peter and Müller (1948): where,   is the Shields number,  c is the critical Shields number, s g is the specific gravity of submerged sediment, and d is a characteristic grain diameter. Here   is calculated as follows: where   2 2 x y U U U , f C is the river-bed resistance coefficient, which is calculated using the Manning-Strickler relation, Here s k is the hydraulic roughness height and h is the flow depth. We assume (Nelson & Seminara, 2012) that the hydraulic roughness height changes linearly with the areal fraction of alluvial cover c p : where sa k is the hydraulic roughness height of alluvial bed material and sb k is the hydraulic roughness height of bedrock itself. The critical Shields number  c is calculated as follows: where  ca is the critical Shields number for a grain on an alluvial bed and  cb is the critical Shields number for a grain on bedrock.
The transverse flux bcn q can be estimated from the relation of Hasegawa (1989), or equivalently that of Mosselman and Crosato (1991):  where   , bs bn u u are the along-streamline and cross-streamline components of the vector of near-bed velocity,    / n is the local riverbed slope in the transverse direction,  s is a coefficient of static friction of sediment (here set to 0.7), and  k is a coefficient of kinematic friction of sediment (0.7). As the flow model calculates the depth-averaged velocities   , x y U U but not the near-bed velocities, the ratio of bs u to bn u is calculated using the following equation, which includes the influence of secondary flow: where  N is a parameter dependent on the strength of the secondary flow and s r is the radius of streamline curvature. Engelund (1974) defines  N as approximately 7.0. Streamline curvature s r is given as follows (see Appendix B): The parameter bs V corresponding to the volume at which the magnitude of bedload transport reaches capacity takes the following form: where s u is the bedload saltation velocity calculated by the following empirical relation (Sklar & Dietrich, 2004)

Bedrock Vertical Erosion
The erosion rate of a bedrock surface is predicted using a 2D version of the model of Chatanantavet and Parker (2009) where  b is the elevation of the bedrock surface and  is the abrasion coefficient. The experimental results of Inoue, Yamaguchi, et al. (2017) showed that the abrasion coefficient for bedrock  can be estimated as follows: where   4 10 C and  T is the rock tensile strength.
The areal fraction of alluvial cover c p is calculated as a function of the sediment layer thickness overlying the bedrock (Inoue et al., 2014;Mishra & Inoue, 2020;Zhang et al., 2015). This method, termed Macro-Roughness Saltation Abrasion Alluviation model (MRSAA) by Zhang et al. (2015), is a descendant of Sklar and Dietrich (2004) which allows routing of alluvium over a bedrock surface.
where mr L is the bedrock macro-roughness (e.g., the standard deviation of bedrock topographic height).

INOUE ET AL.
10.1029/2020JF005645 7 of 19 Bed surface elevation variation is given as follows: where  is the elevation of the top of the river bed, and the subscripts a and b represent alluvium and bedrock, respectively. We assume that the debris of bedrock abrasion is small enough to be treated as wash load, and do not consider an increase in sediment supply caused by bedrock bed and bank abrasion in our model.

Bedrock Lateral Erosion
Various studies suggest that bedrock bed erosion by saltating bedload depends purely on the rate of particle collision with the bedrock (Inoue et al., 2014; J. P. L. Johnson & Whipple, 2010;Sklar & Dietrich, 2004).
Hence the erosion rate of the bedrock bed strongly depends on the magnitude of bedload transport rate. Conversely, studies have suggested that bedrock bank erosion may not share a similar dependency on the magnitude of bedload transport rate. For example, in a naive linear channel model, where the bedload only travels streamwise, the bedload saltates solely parallel to the sidewalls, thus reducing the rate of sidewall collision to zero. Also, the experimental results of Mishra et al. (2018) have shown that the lateral erosion rate is almost zero in a straight channel with a bed fully covered by uniform sediment. They suggested that lateral gravel movement caused by flow meandering and bar formation controls the lateral erosion of bedrock walls. As previously mentioned, the lateral bedload transport rate is modeled as zero at the bank in an Eulerian dispersion model because the number of colliding particles and the number of rebounding particles are equal. However, in reality, laterally moving particles will still impact the wall, causing erosion. Therefore, in this study, we propose a simplified version of the model proposed by Inoue, Parker, et al. (2017) that relates the lateral bedrock erosion rate to the lateral sediment transport rate measured at a distance b L away from the wall; where,  bank is the sidewall bedrock abrasion coefficient with the same dimensions as the abrasion coefficient for bedrock bed  (m −1 ), and b L is the thickness of a near-bank boundary layer over which the lateral bedload transport rate drops to zero. Here, we use the term boundary layer in a general way to characterize sediment transport transitioning from a relatively large value at a distance b L from the sidewall to zero at the sidewall boundary. We do not refer to a turbulent boundary layer that develops in flows in proximity to a fixed boundary such as the bank. In this study, we assumed that  bank is equal to . The thickness of the boundary layer b L may change depending on hydraulic conditions, but it has not been clarified to date. For simplicity in these simulations, we set b L as 10% of the channel width.
In natural rivers, however, irregular movement of particles is driven by the effect of bed roughness (Fuller et al., 2016). Here we include the roughness effect in a simple way by introducing a specified coefficient; where  r characterizes the probability of collision with the sidewalls due to irregular movement of particles transported in the longitudinal direction in the mean. the heuristic for randomness in particle trajectories driven by bedrock topographic roughness. However, since accurate prediction of hydraulic roughness should take into account not only the bedrock topographic roughness but also the arrangement of bed unevenness (Mishra & Inoue, 2020), we do not link the hydraulic roughness height s k and the roughness effect coefficient  r in this study.
In our model, sediments move toward left bank when by q is positive, and sediments move toward right bank when by q is negative. Equation 21 needs to be divided into two cases each for the left and right banks:  10.1029/2020JF005645 9 of 19 Figure 10. Time variation in Run 2 of (a) channel sinuosity, (b) mean bed slope, (c) mean channel width, (d) bedload transport rate near the right bank at high flow averaged along the streamwise direction, and (e) lateral sediment transport rate toward the right bank at high flow averaged along the streamwise direction.

Land Accretion
In alluvial rivers, as bar height increases, vegetation grows and traps suspended sediment and the inner part of a bar eventually is converted from river bottom to adjacent land. Submergence then becomes very rare, occurring only during large flood events. In order to include this process in numerical calculations, Asahi et al. (2013) removed the dry area from the computational grid and redistributed the calculation area. Even in mixed alluvial-bedrock rivers, vegetation grows on bars composed of gravel. In addition, in mixed alluvial-bedrock rivers, outside bedrock benches are left behind when the channel incises downward. Inoue, Parker, et al. (2017) removed bedrock benches from the computational grid where the angle of the bench-thalweg juncture exceeded a specified angle, but this method cannot remove a point bar when the lateral slope is relatively gentle. Therefore, in this study, we remove the area where the water depth is below a specified depth dry h from the computational grid (Figure 3), and redistribute the calculation area as per the method proposed by Asahi et al. (2013).
In the model presented here, the land outside the channel is not modeled, and the flow does not go over the sidewalls. As previously mentioned, we assume that the elevation of the area that accretes land increases due to trapping of suspended sediment during low flow, and is then no longer inundated even during high flow. The land returns to the river when the sidewall is laterally eroded.

Calculation Conditions
The planform shape of the initial channel is set using a sine-generated curve. In all simulations, the initial channel width is 20 m with 10 grid cells, initial bed slope is 0.01, initial meandering wavelength is 450 m with 24 grid cells, initial meandering amplitude is 50 m, initial sinuosity is 1.03, and grain size is 20 mm. The initial bedrock surface elevation is uniform in the transverse section, and the initial alluvial layer is spread uniformly across the bedrock surface. As we use cyclic boundary conditions, such that, water and sediment leaving the downstream end of the channel directly enters the upstream end. The total volume of alluvium placed in the initial channel is seen to be a conserved quantity. In our model, sediment supply and flow discharge are not completely decoupled because the sediment transport rate (and also the sediment supply rate) is calculated using the fluid shear stress, which changes with flow discharge. When the initial sediment volume is larger than the transport capacity, the sediment supply rate is equal to the sediment transport capacity. When the initial sediment volume is less than the transport capacity, the sediment supply rate is limited. This situation is generally observed in natural bedrock rivers. Being able to calculate supply-limited conditions is a strength of our model. Under cyclic boundary conditions, the variation of bed height and channel width is the same at the upstream end and the downstream end because the water depth, flow velocity, and sediment transport rate are forced to be the same due to cyclic boundary conditions. Base level, on the other hand, is not controlled in the version of the model presented here with cyclic boundary conditions. In general, the hydraulic roughness height of an alluvial bed is taken to be  1 4d. Here, the alluvial hydraulic roughness sa k is set to 2.5d (median of 1d to 4d). The value of the critical Shields number for an alluvial bed is generally 0.03-0.073 (Buffington & Montgomery, 1997). This model uses 0.05 for  ca . To simplify our simulations here, the hydraulic roughness height for bedrock sb k and critical Shields number for grains over bedrock  cb are taken to be equal to those of an alluvial bed sa k and  ca , respectively. The bedrock macro-roughness length mr L is set as 1 m. The rock tensile strength  T is taken to be 0. based on measured values of weak bedrock outcrops in the Shikaribetsu river and the Kushiro River (Hokkaido, Japan) by Inoue, Yamaguchi, et al. (2017).
We change the initial alluvial thickness and the flow discharge in order to explore the sensitivity of the model (Table 1). In order to evaluate the distribution of the lateral sediment transport rate when alluvial point bars reach a dynamic equilibrium state, the alluvial layer is deformed in Run 1 but erosion was set to zero in both the bedrock bed and walls. In Runs 2-11, we consider bedrock erosion of both bed and walls, and investigate the morphodynamics of meandering incised into bedrock. In Runs 2-11, we use cyclically varying hydrographs, because land accretion of alluvial point bars mainly occurs at low flow.
Understanding long-term river evolution demands a computationally expensive or unwieldy consideration of varying discharges and their sequences over a time scale of at least centuries. Thus, in this study, we use a simple repetitive two-stage hydrograph ( Figure 4) where high discharge emulates large flood events and low discharge promotes land accretion. The low discharge time period is substantially cut short so as to calculate only the time period over which land accretion occurs, whereas the high discharge time period is calculated as is. The total calculation time is 7200 h for all runs (240 cycles of two-stage hydrograph). Assuming that one combined cycle of high flow and low flow corresponds to 1 year, the calculation time amounts to 240 years.
In Runs 2-6, the initial alluvial thickness is specified in terms of five values (0.5, 1.0, 2.0, 5.0, and 10.0 m). As the boundary conditions are cyclic, this thickness serves to specify the sediment supply, which then remains constant throughout a numerical run. In Run 2 and 7-10, five values (50, 100, 150, 225, and 300 m 3 /s) are assigned to the high flow discharge of the repetitive two-stage hydrograph. In Runs 2 and 11, we investigate the effect of collision probability of irregular particle motion caused by the bed roughness during channel evolution.
The critical flow depth for land accretion dry h is set as 1% of the average flow depth at the low flow discharge for Runs 2-11 ( dry h = 0.05 m). Figure 5 shows an example of land accretion during Run 2. The high flow depth before land accretion is roughly 0.33-0.66 m (Figure 5a). We assume that the elevation of this area increases due to the deposition of suspended load during low flow (Figure 3b), and the area is not inundated significantly even during high flow (Figure 3c).

Lateral Sediment Transport Rate
In Run 1, for which bedrock bed and bank erosion is not permitted, the height of alluvial point bars reaches an equilibrium state after roughly 4 h ( Figure 6). Point bars form on the downstream side of the meander bend apex, and the river bed elevation tilts downward from inner bank to outer (Figure 7a). This transverse bed slope should force the gravel toward the outer bank, but the secondary flow near the river bed suppresses this movement of gravel (i.e., Equations 11 and 12). In case of a uniform bend, that is, one with constant curvature, Inoue, Parker, et al. (2017) show that channel cross-section morphology adjusts to an equilibrium state so that the effects of the lateral bed slope and the secondary flow cancel each other, and the lateral bedload transport rate becomes zero unless lateral erosion is initiated by an increased sediment supply (e.g., Figure 12 in Inoue, Parker, et al., 2017). In a sine-generated channel, however, the lateral bedload transport does not become zero even if the alluvial point bars reach equilibrium state (Figure 7b) because there is a phase lag between riverbank lines and streamlines caused by fluid inertia (Figure 7c). In this condition, the streamline-perpendicular bedload transport rate bn q is zero, but the bedload transport rate perpendicular to the river bank line (or computational grid line) by q is not zero, because by q includes the bedload transport rate in the streamline direction bs q as shown in Equation 5. This means that as long as there is a phase lag between the streamline and the riverbank line, gravel collides with bedrock sidewalls and erodes them.

Evolution of Bedrock Incised Meanders
In Run 2 with a cyclic hydrograph (repeated 150-15 m 3 /s), the channel cuts downwards and outwards while meandering (Figures 8 and 9). Lateral incision is computed in terms of the intensity with which streamlines converge toward a bank, driving the grains to strike it. Although the increase rate in sinuosity decreases gradually, the channel sinuosity does not become constant (Figure 10a). Decreasing the bed slope and increasing the river width reduces the magnitude of sediment transport (i.e., tools for bedrock erosion), and consequently decreases the increase rate in channel sinuosity (Figures 10b-10d). On the other hand, the following two factors drive the lateral bedrock erosion and inhibit arrival of the equilibrium state: (1) the lateral sediment transport rate decreases with the decrease in the magnitude of sediment transport rate but does not vanish (Figures 10e and 11) because the phase lag between the streamline and the bank line is not eliminated even if the meandering progresses; (2) fluctuations in water depth change the strength of the secondary flow and repeatedly increase lateral sediment transport toward the sidewalls, in particular, during recession limbs in which the secondary flow is weak ( Figure 11). Thus, the sinuosity does not reach constant, and the meander amplitude slowly increases while the channel incises downward, as shown in Figure 9.

Effects of Sediment Supply, Flow Discharge, and Bed Roughness
Comparing the results of Runs 2, 3, 4, 5, and 6, we can obtain a picture of the effect of sediment supply on bed evolution (Figure 12a). We can see that the vertical erosion depth decreases with increasing sediment supply, because the area of bedrock exposure decreases. When the initial alluvial thickness exceeds 2 m, the bedrock is almost covered with gravel, and the vertical bedrock erosion is almost zero. On the other hand, the channel width and sinuosity increase with increasing sediment supply (Figures 12b and 12c). When the sediment supply is relatively low (e.g., Run 3), bedrock benches form at the outer bank (Figure 13b) because the width of the point bar is smaller than the channel width (Figure 13a). These benches turn into land as the thalweg incises downward, after which the banks of the bench-thalweg are eroded laterally ( Figure 13c); That is, the formation of the bench reduces channel width and stagnates the growth of sinuosity. When the sediment supply is very large (e.g., Run 6) the bedrock is fully covered by alluvium, no bedrock erosion occurs in the vertical direction, and only the uncovered sidewalls are eroded laterally (Figure 13d).
Similarly, by comparing the simulations of Runs 2, 7, 8, 9, and 10, we can study the effect of variation in the high flow discharge on channel evolution (Figure 14). It appears that vertical erosion increases with increasing flow discharge, because high discharge leads to a greater transport capacity, in turn leading Changes in bedrock elevation overlain every 600 h in Run 3; (d) Changes in bedrock elevation overlaid every 600 h in Run 6. A lower sediment supply leads to larger vertical erosion, as well as river width reduction with bench formation. In the case of a high sediment supply, the channel migrates laterally without vertical erosion.
to reduced alluvial cover. Also, high discharge contributes to significant lateral sediment transport, resulting in high channel sinuosity. Compared to the above two relationships, the relationship between channel width and flow discharge is more complex. More specifically, the channel width is maximum at a specific flow discharge. While ever higher flow causes ever larger lateral erosion, it also increases the height of the point bars, subsequently promoting land accretion at low flow ( Figure 15). Focusing on the shape of the meanders, when the high flow discharge is relatively low, a meandering shape with upstream skewness appears (Figure 16).
When stochastic particle collisions caused by riverbed roughness are accounted for, the channel sinuosity increases by 24% and vertical erosion increases by 19% (Figure 17). In addition, the channel shape shows clear upstream skewing, as is commonly observed in well-developed alluvial meanders (Guo et al., 2019).

Discussion
In a meandering bedrock river, even when the river bed is in a dynamic equilibrium state, gravel can be moved laterally toward bedrock sidewalls due to the difference in direction between the streamlines and the bank lines, so resulting in sidewall erosion. More specifically, sidewall bedrock erosion occurs where the bedload collides with the sidewall (Figures 7  and 8). The lateral erosion rate depends on sediment supply and flow discharge, and channel sinuosity increases with an increase in these parameters (Figures 12c and 14c). In particular, an increase in sediment supply simultaneously suppresses vertical erosion, but does not suppress lateral erosion (Figures 12a and 13d). Vertical erosion is dominant with small sediment supply (Figures 12a and 13c). These results suggest that variation of sediment discharge can result in the formation of strath terraces, an observation which is macroscopically consistent with the findings of previous studies (Fuller et al., 2009;Gilbert, 1877;Hancock & Anderson, 2002). That is, changes in sediment and flow discharges induced by, for example, climate change will cause a change in the ratio of vertical to lateral erosion of bedrock rivers.
Looking at the calculation results in detail, however, the ratio of vertical to lateral erosion rate is seen to vary in time even if the sediment and flow discharges are constant. In Figure 9, the ratio of vertical to lateral erosion rate near the bend apex is roughly 0.04 from 0 to 600 h, and roughly 0.3 from 6600 to 7200 h. Increasing channel bend amplitude strengthens the secondary flow and reduces the rate of particle collision against the bedrock sidewalls (see Figure 10e). A decrease in bed slope (see Figure 10b) decreases the sediment transport rate (see Figure 10d), so reducing the availability of tools for bedrock bank erosion. This suggests that sinuosity may eventually stabilize (see Figure 10a) as the tools decrease, especially when there is a threshold strength for bedrock erosion. In addition, although decrease in sediment supply decreases lateral erosion and increases vertical erosion, changes in channel width have a role in suppressing the influence of sediment supply on lateral erosion. Channel narrowing promotes lateral erosion by increasing contact of sediment with the sidewalls of the point bars and, at the same time, suppressing vertical erosion due to increased alluvial cover in the narrower channel (Figure 12a). Channel widening with large sediment supply has the opposite effect. These results indicate that in addition to sediment supply, channel morphodynamics, such as changes in channel bend shape, bed slope, and channel width are also important factors controlling the ratio of vertical and lateral erosion. The above findings support the argument of Finnegan and Dietrich (2011) that terraces may form without changes in boundary conditions. Bedrock meanders with clear upstream skewness form in Run 7 (Figure 16a) and Run 11 (Figure 17b). In Run 7, the effects of land accretion are larger than those of bank erosion because of the relatively small flow discharge. In Run 11, bank erosion includes the effects of stochastic collisions of gravel transported along bank with the sidewall. Compared to the lateral sediment transport rate, the longitudinal sediment transport rate is simply dependent on the fluid shear stress, which is a major factor in alluvial bank erosion. Guo et al. (2019) indicate that upstream skewness tends to be observed in alluvial channels with high sinuosity. However, our modeling suggests that a bedrock channel may not show upstream skewness when the effects of land accretion and stochastic collisions are relatively small.
In this study, because a gentle meander is given as an initial condition and a periodic boundary is used, the meander shapes evolve regularly, except in Run 3 with low sediment supply. In the regularly meandering channels, the channel is wide on the upstream side of the bend apex and narrow on the downstream side (e.g., Figure 18a) because alluvial point bars are more well developed on the downstream side. Similar point bar development has been observed in the experiment performed by Mishra et al. (2018). In addition, the bed slope is steep near the bend apex ( Figure 18b). In contrast, in the irregularly meandering channel, the channel width does not change significantly ( Figure 18c) and the bed slope is irregular along the longitudinal direction (Figure 18d). In Run 3 with a low sediment supply, bedrock benches form near the outer banks and complicate the flow fields. We expect this to be a factor related to irregular meander shape.
Cutoffs are not modeled here, and indeed, our calculations are not performed long enough for adjacent channels to overlap in this study. We believe that cutoffs can occur if irregular meandering like Run 3 is allowed to evolve further. Since the lateral erosion rate is low in Run 3 with low sediment supply, it will take a long time for a cutoff to occur. Therefore, in order for a cutoff to occur, two "eras" may be required: an "era" with low sediment supply in which the meandering shape becomes irregular, and an "era" with high sediment supply in which the meandering amplitude strongly increases. The relationship between meander shape and temporal change in sediment supply is an appropriate subject for future research.
When a bedrock side-wall overhangs and then collapses, the bedrock will deposit as large blocks on the riverbed. The bedrock blocks can "cover" the bank and inhibit bank erosion, like slump blocks of cohesive material in an alluvial river (Parker et al., 2011). However, the bedrock blocks can also be "tools" that promote downstream bank erosion after the bedrock blocks weather and turn into coarse sediment. The "cover" and "tools" effects of bedrock slump blocks are not included in our model, and instead are a topic for future modeling.
Other future research topic is whether meandering develops if we start with a straight channel. Recent studies show that alternate bars or longitudinal grooves form on a straight bedrock channel (Chatanantavet & Parker, 2008;Inoue et al., 2016;Inoue & Nelson, 2020). These bed forms may deflect the flow and trigger the development of meandering. As there is still debate about whether alternate bars can trigger meandering in alluvial channels (e.g., Seminara, 2006;Whiting & Dietrich, 1993, with particular emphasis on the role of bar-bend resonance), this represents very difficult and exciting challenge in the future study of bedrock channels.

Conclusions
It has long been known that mixed bedrock-alluvial rivers, that is, bedrock rivers with a partial alluvial cover, can show as intense a pattern of meandering as purely alluvial rivers. The study of the bedrock-alluvial case has lagged behind the alluvial case, largely due to the difficulty of field observations. Here we show the results of a numerical model that accounts for abrasion-driven meandering processes in a physically realistic way, as follows.
• The total amount of alluvium is conserved, using cyclic boundary conditions. • Vertical incision is computed in terms of abrasion of the bedrock by saltating grains. • Lateral incision is computed in terms of the intensity with which streamlines converge toward a bank, driving the grains to strike it. Allowance is also made for fluctuations in the paths of grains moving parallel to the bank in the mean. • A two-step hydrograph accounts for the role of flow variation. Land accretes where the depth falls below a specified value at low flow.
The model is able to replicate the evolution of high-amplitude meandering over a wide range of conditions. The rate of sinuosity growth declines in time, but the bends do not stabilize. For the same hydrograph, increasing the initial mean thickness of alluvium suppresses vertical incision and increases the growth rate of sinuosity. At a sufficiently high thickness, the channel migrates laterally without bed lowering. Bend shape can change dramatically with increasing alluvial thickness, with high thickness favoring more regular bend trains. For the same initial mean thickness of alluvium, on the other hand, increasing the peak discharge increases the rate of vertical incision and rate of growth of sinuosity.

INOUE ET AL.
10.1029/2020JF005645 15 of 19 Figure 17. Changes in bedrock elevation overlain every 600 h in Run 2 (without irregular particle motion) and Run 11 (with irregular particle motion). In Run 11%, 10% of the sediment transport rate in the longitudinal direction collides with the sidewalls due to irregular movement caused by bed roughness. Upstream bend skewing clearly appears only when the influence of factors that control alluvial meandering (e.g., land accretion, shear stress) is large. Most of the modeled cases do not show clear upstream skewing, a feature that replicates the planform of the Shimanto River shown in Figure 1b.
The model invites relatively straightforward extension to problems of long-term geomorphic significance. One such problem is the formation of strath terraces. These can form, for example, due to long-term changes in sediment supply. Although computationally expensive, the model presented here can handle such variation in a straightforward way. where f C is the river-bed resistance coefficient.
The parameters x D and y D expand into where  t is a coefficient of eddy viscosity given by where  is the von Kármán constant and  u is the shear velocity. Figure B1. The relationship between  s n coordinate system and  x y coordinate system.