Steady‐State Parallel Retreat Migration in River Bends With Noncohesive (Composite) Banks

A substantial body of research has addressed the equilibrium cross‐sectional geometry of straight noncohesive channels, along with bends having fixed outer banks. However, development of a characteristic cross‐section during active migration has been confounded by inaccurate treatment of noncohesive bank erosion processes. This analysis characterizes a steady‐state migrating cross‐section and the associated migration rate for the highly conceptualized case of an infinite bend of constant centerline radius with noncohesive lower banks consisting of uniform‐sized grains mobilized as bedload. Analytical, numerical, and field analyses are presented to rationally constrain the geometry and obtain a physically based migration rate equation dependent on the following dimensionless groupings: excess Shields stress, flow depth to radius of curvature ratio, and noncohesive bank thickness to grain size ratio. Migration rate is shown to be dictated by transverse sediment flux at the thalweg due to secondary flow, not bank slope as in previous formulations developed from similar principles. Simple outward translation of the noncohesive portion of the cross section is achieved with concurrent mass failure in the upper cohesive layer. Such translation cannot be achieved using a linear excess shear formulation applied to noncohesive fluvial erosion. A numerical model of cross‐sectional evolution to steady‐state migration is developed; when applied to the lower Mackinaw River in Illinois, it reveals that the river behaves as if the critical shear stress is considerably larger than that indicated by the grain size distribution. This conceptualized treatment is intended to provide a canonical basis of comparison for actual meander bend geometries.

. In the past several decades, these processes have been incorporated into a multitude of bank erosion models of varying sophistication; descriptions of relevant physical processes and many of the modeling approaches can be found in  and Klavon et al. (2017).
Properly characterizing the fluvial erosion of noncohesive bank materials remains a weakness of morphodynamics models. This weakness exists despite substantial attention being given to the topic of noncohesive bank deformation during the 1980s and 1990s (e.g., Diplas, 1990;Ikeda, 1981;Ikeda et al., 1988;Kovacs & Parker, 1994;Pizzuto, 1990; Thorne et al., 1998). These studies characterized evolution of banks to a threshold channel configuration when a straight channel is subjected to excess boundary shear stress with the bank material mobilized as bedload. The threshold channel concept (Glover & Florey, 1951), as modified by Parker (1978) to account for an active bedload transporting region, provided a coherent framework upon which bank deformation results could be generalized. Those studies that incorporated numerical modeling implemented several approaches to account for momentum extraction due to the transverse boundary layer associated with the bank; the resulting boundary shear stress distribution was then used to drive the boundary deformation due to transverse bedload flux divergence terms in the sediment mass conservation equation (Exner equation). Despite the success of these early modeling approaches, the numerical techniques developed for straight channels has been slow to be implemented into morphodynamics models for natural channels. This can be attributed primarily to the increased complexity of natural channel alignments and the consequent increased computational expense of accurately characterizing boundary shear stresses at a suitable spatial resolution, an issue highlighted by Rinaldi and Nardi (2013).
Presently, the most commonly used reach-scale modeling techniques that address fluvial erosion of basal noncohesive bank materials can generally be classified into two broad categories, which depend on the bank erosion modeling framework: (a) models that treat bank erosion through analysis of individual cross sections, in which the boundary shape and boundary shear stress distribution can be discretized at a denser transverse resolution than the width of the entire bank; and (b) spatially continuous two-dimensional (2D) morphodynamics models, in which the entire width of the bank at any longitudinal position along the stream may be comprised by few (or one) numerical cells.
With respect to the first category of bank erosion models, the most common method for fluvial erosion specifies the local fluvial erosion rate in the direction normal to the surface, E, as a linear function of the excess boundary shear stress: where k is an erodibility coefficient, τ b is the boundary shear stress vector, and τ c is the critical boundary shear stress. This linear formulation is also known as the Ariathurai-Partheniades equation (Ariathurai, 1974;McAnally & Mehta, 2001). The linear excess shear stress method is used in models, such as BSTEM (Simon et al., 2000), CONCEPTS (Langendoen & Simon, 2008), and various physics-based models developed for academic research (e.g., Darby et al., 2007;. However, the formulation was developed based on experimental results using cohesive sediments. A recognized knowledge gap has been identified regarding suitable ways to specify the erodibility coefficient for noncohesive bank materials (e.g., Rinaldi & Nardi, 2013).  were compelled to assume that the equation is robust, given the lack of a suitable alternative, but they acknowledge that the physical basis of the equation is questionable. Waterman and García (2016) described how Equation 1 is an entrainment rate equation that applies to sediments transported at less than the transport capacity of the flow, which improperly characterizes sediment mass conservation when dominated by noncohesive sediment transport. Nevertheless, the excess shear stress formulation continues to be extensively used in the cross-sectional models, primarily due to its ease of numerical implementation and long history of usage.
With respect to the second category of models, in which the bank region is discretized with few numerical cells, simplified conceptualizations of fluvial erosion and mass failure are generally "lumped" into algorithms of varying physical realism (Stecca et al., 2017). This lumped (or bulk) treatment is a practical necessity due to computational limitations with respect to modeling reach-scale spatially continuous domains. With this limitation in mind, the accurate expression of sub-grid scale bank processes as net bank retreat remains one of the major challenges of modern morphodynamics modeling (Siviglia & Crosato, 2016). The simplest of these lumped bank erosion models, which are not specific to either cohesive or noncohesive banks, specify the bulk bank retreat rate as a linear function of either the excess near-bank velocity or the excess near-bank boundary shear stress.

Motivation, Objectives, and Approach
The bank-integrated Exner equation approach for bank migration requires calculation of the transverse flux of sediment from the base of the bank region, which is dependent on both the effects of secondary flow and gravity due to bank slope on bedload movement. The integral treatment as implemented by Hasegawa (1989) and Parker et al. (2011) requires a similarity function for the cross-sectional shape of the bank. In those analyses, the bank shape was assumed to have a uniform transverse slope; however, a rational, physically based manner of specifying the bank shape a priori has never been established. A characteristic cross-sectional shape for a meander bend does not exist in the same sense as for an idealized straight channel. The width frequently changes due to differential rates of outer bank erosion and inner bank accretion (e.g., Eke et al., 2014;Mason & Mohrig, 2019;Zhao et al., 2021). Slump block armoring of the bank alters bank shape and erosion rates (e.g., Hackney et al., 2015;Motta et al., 2014;Parker et al., 2011). The bank shape evolves depending on the discharge (e.g., Zhao et al., 2020), which is variable, and the velocity field also varies through the bend. Thus, the idealization necessary to ascertain a characteristic shape will never strictly represent a real meandering river. Nevertheless, this analysis is motivated by the primary goal to identify a canonical case of meandering river cross-sectional shape to inform the shape constraint required in the bank-integrated Exner equation approach. Whether the results from the idealized case are suitable for practical use ultimately depends on the degree to which they correspond to real rivers; but even in those cases where they do not, they will provide a basis for comparison to highlight similarities and differences to shed light on those physical factors responsible for the differences.
The following simplifications are implemented in the current analysis: (a) fully developed flow in a constant radius bend of low to moderate curvature based on the nomenclature of Blanckaert (2011); (b) uniform-sized sediment mobilized only as bedload; and (c) the fine-grained upper layer plays no role in the outer bank migration dynamics through armoring of the noncohesive layer by cohesive material that undergoes mass failure. The cross-sectional shape is sought that migrates without changing shape (parallel retreat) under a constant discharge. More physically realistic representations of all the complexities of natural channels can be developed in the future once the simplest condition is better understood. Note that the assumptions of fully developed bend flow and a relatively time invariant cross-section does not represent transitional and braided rivers. Such rivers still need to be treated differently than meandering rivers with composite banks in numerical models.
With these simplifications, the noncohesive layer is effectively treated independently from the cohesive top layer, which warrants justification. The mass failure of the cohesive upper layer and its armoring effect is clearly a dominant factor in determining the migration rate of rivers with composite banks. However, the influence of 4 of 33 slump blocks will be dependent on the thicknesses of the upper cohesive layer relative to the underlying noncohesive layer. Parker et al. (2011) referenced findings from Muddy Creek, Wyoming (Dietrich & Smith, 1983) to illustrate the importance of slump block armoring. Muddy Creek had a bankfull depth of ∼0.5 m, and failed slump blocks were large enough to protect the entire lower bank. On the other hand,  and  describe a bend of the Wabash River in which the lower 11 m of the bank region consisted of sand/gravel, while the upper cohesive layer was generally less than 2 m thick. In such a case, slump blocks would be expected to protect only a small portion of the bank region. Therefore, the influence of cohesive slump blocks can be bounded by two limiting cases: (a) the case where slump blocks have no influence on the migration rate; and (b) the case where slump blocks cover the entire bank and erode at rates orders of magnitude slower than the noncohesive layer, and whose resistance to erosion entirely controls the migration rate. The idealized case considered describes one extreme that bounds the range of possibilities, and deviations from the idealized case will shed light on the influence of slump blocks on channel shape and migration rate.
The first objective of the analysis is to suitably constrain the bank shape for use in the integral treatment. The second objective is to utilize the shape constraint to express the bank-integrated Exner equation for migration rate in meandering rivers with composite banks. The third objective is to demonstrate numerically that parallel retreat as conceived in the analytical treatment as a simple outward translation of the noncohesive portion of the cross-section (with concurrent mass failure in the cohesive upper layer) is feasible under the simplifications considered.
The paper is organized as follows. In Section 2, the integrated Exner equation is explored to identify the conditions required to achieve parallel retreat and to establish the basic form of the migration rate equation whose terms must be quantified. Section 3 theoretically establishes the parallel retreat shape constraint of zero transverse slope (horizontal) at the channel thalweg; this is contrasted with the implied similarity functions of Hasegawa (1989) and Parker et al. (2011). The new shape constraint allows further quantification of terms in the migration rate equation, which is then expressed in dimensionless form as a function of the excess Shields stress. In Section 4, a simplified numerical model is presented to demonstrate that parallel retreat as simple translation of the noncohesive portion of the cross-section is feasible as conceived in the theoretical analysis. In Section 5, the numerical model is applied to a bend of the Mackinaw River, and the simulation parameters that yield results that best represent the field observations are described. Finally, in Section 6 organizing principles are discussed that provide a coherent framework upon which noncohesive bank deformation can be generalized across both straight and meandering reaches; the results are also discussed in terms of relating the highly simplified conceptual model to the complexities of real rivers.

Integrated Exner Equation and Assumptions
The infinite bend configuration (fully developed flow, constant centerline radius) used throughout this analysis is illustrated in Figure 1. The governing equations are most conveniently expressed in modified cylindrical coordinates (s, r, z), where s is the streamwise coordinate, r is the radial coordinate, and z is the vertical coordinate. The modified cylindrical coordinates are converted to standard cylindrical coordinates based on the relationship s = rθ, where θ is the angle turned about the center point in a plane where z = constant. The longitudinal slope of the channel is considered small, such that (s, r, z) are orthogonal coordinates. The variables r in , r c , and r out represent the radii of the inner boundary, the centerline, and the outer boundary of the channel, respectively. The radial position of the center point of the circle is denoted r cp , which always equals 0.
The constant migration rate M is sought for a channel geometry that does not change during migration under a steady bankfull discharge Q bf . To maintain constant radial values as the channel migrates, the (s, r, z) coordinate system must represent a translating frame of reference rather than concentric expansion on a fixed frame of reference. The resulting governing equations have the same basic structure as the more commonly used intrinsic (s, n, z) coordinate system (e.g., Parker et al., 2011) where the origin of the transverse n coordinate is associated with the dynamic channel centerline. In Figure 1 inset, the translating frame of reference is illustrated with respect to a fixed frame of reference, in which the radial coordinate for the latter is denoted with R. For each value r cp , r in , r c , and r out in the translating system, an associated value R cp , R in , R c , and R out , respectively, exists in the fixed frame of reference. The following relations describe the time-invariant radial quantities under the constant migration rate M.
The constant channel width B is expressed as B = r out − r in . The steady discharge Q bf is entirely conveyed between r in and r out . The volumetric bedload transport rate per unit width vector, q, is decomposed into a streamwise component, q s , and a transverse component tangent to the boundary surface in the r-z plane, q rt . In modified cylindrical coordinates, the Exner equation is expressed as: where η represents the surface elevation above an arbitrary vertical datum and λ is the porosity of the granular medium. The partial derivative of q rt is taken with respect to the coordinate r rather than a boundary-fitted coordinate defined tangent to the surface in the r-z plane, as derived in Waterman and García (2016). In steady, uniform bend flow, the divergence terms in the s-direction vanish, allowing simplification to: Utilizing the definition of the total derivative of η(t,R) yields the following: For parallel migration, in which the shape does not change within the translating frame of reference, dη/dt = 0. This treatment presupposes that the parallel retreat process is not a time-averaged depiction of a cyclical process of basal erosion and mass failure in the bank region, but rather a simple translation of the entire cross section due to local erosion and deposition. Substituting from Equation 4, and noting that ∂R/∂r = 1 from Equation 5, Equation 8 is re-expressed accordingly: = − Equation 9 ensures that the bank region is erosional when considered on a fixed frame of reference due to ∂η/∂r being positive. Substituting Equation 9 into Equation 7 yields: Following the technique of Hasegawa (1989) and Parker et al. (2011), Equation 10 is integrated in the radial direction; the limits of integration can be set at any arbitrary r positions on the cross-section. Integrating over the r-direction between arbitrary radial positions r A and r B yields: Note that the same result is obtained by integrating Equation 7 over the r direction and implementing the Leibniz integral rule, which was the technique used by Parker et al. (2011). Ultimately, Equation 12 will be used to formulate the migration rate; a key aspect of the analysis is to appropriately establish the positions r A and r B where q rt and η can be readily quantified.

Boundary Conditions and Their Implications
The boundary condition is imposed that q rt | rout = 0; that is, no bedload enters the domain from the cohesive upper layer as the channel migrates. Under the idealized condition, the upper cohesive layer passively fails as the noncohesive layer migrates, and that material is carried away in suspension. The boundary condition is also imposed that q rt | rin = 0. The basis of the latter boundary condition is that, by problem definition, no discharge is conveyed where r < r in . Thus, the mean flow velocity, bed shear stress, and bedload transport rate at the boundary are equal to zero. The conceptual interface at r = r in is illustrated in Figure 2; the flow depth is non-zero at the interface, but dense vegetation in the region r < r in is presumed to maintain negligible velocity. The necessity of the non-zero flow depth at the inner boundary is described later in this section. Figure 2 shows that η| rout is the top elevation of the noncohesive layer at the interface with the overlying cohesive layer, and r thal is the radial coordinate of the thalweg. Figure 2 illustrates regions of aggradation (net deposition) and degradation (net erosion) that must exist in a fixed spatial frame of reference under conditions of parallel retreat. The shape of the cross-section in Figure 2 is not imposed; the only shape assumption made is that ∂η/∂r ≤ 0 throughout the region left of the thalweg and ∂η/∂r ≥ 0 throughout the region right of the thalweg. The cross-sectional shape illustrated is only provided to aid the description of the problem setup.
With the shape assumption that ∂η/∂r ≤ 0 throughout the region left of the thalweg, the outermost location of maximum η must be aggradational and emplaced as the channel migrates. Since aggradation necessitates the condition that ∂q rt /∂r ≠ 0, the flow must be competent to transport bedload infinitesimally close to the location of maximum η. This establishes the non-zero flow depth requirement at the inner boundary illustrated in Figure 2.
The q rt | rout = q rt | rin = 0 boundary conditions also have implications on the elevation of η| rin relative to η| rout . To achieve parallel migration, the volume of noncohesive sediment eroded must be exactly balanced by the volume deposited in the fixed frame of reference. Dividing the radial slice ( Figure 1) into numerous control volumes in the radial direction, the radial boundaries have increasing arc lengths toward the outside of the bend. Thus, the eroded sediment originates from a region with larger control volumes than the region where sediment is deposited. To achieve the volume balance, deposited sediment must reach a somewhat higher elevation at the inner edge than the elevation of eroded material at the outer edge to compensate for the control volume shape. The mathematical demonstration of this assertion is provided in Appendix A. This requires that η| rout cannot be arbitrarily high on the outer bank but must be lower than the elevation at which bedload can be transported at the inner bank.

Development of the Shape Constraint
The key to using Equation 12 is appropriately establishing the positions r A and r B where q rt and η can be readily quantified. At the outer boundary, q rt = 0, and η| rout can be estimated with respect to the top of bank using a typical depth of the cohesive top layer. Therefore, r out is the most straight-forward choice for setting r B in Equation 12. Selecting the position r A is more challenging, as q rt is a function of the covarying hydrodynamics and boundary shape, which are not known a priori. A rational shape constraint to generalize the expression of q rt at a specific location r A is therefore sought. The following subsections describes the formulation of q rt on a transverse slope in the presence of secondary flow and describe the importance of properly establishing the transverse slope at the The spatial variables identified on the surface are at time t; η represents the surface elevation. The radial component of bedload transport rate tangent to the surface (q rt ) is also illustrated. The circled numbers indicate the following characteristics of the regions that form the channel boundary and just beyond it: (1) cohesive upper layer; (2) noncohesive lower layer; and (3) vegetated over-bar area that is submerged under moderate and high flow conditions. inside limit of integration (r A ). Next, the thalweg position is deemed the most logical position to set r A , and the transverse slope of that local minimum of η is explored. The thalweg is shown to have a transverse slope equal to zero, which allows a straightforward expression of q rt at that position for use in Equation 12. The presence of an abrupt slope discontinuity at the thalweg is shown to be unlikely under the assumptions of this analysis. Figure 3 is a definition diagram showing variables required for the bedload transport rate evaluation, and previously undefined variables shown are described after the figure.

Formulation Used for q rt in the Presence of Secondary Flow
On Figure 3, the variable δ is the deviation angle of τ b relative to the s-direction, with the positive direction always toward the outer bank; β is the deviation angle of q relative to the s-direction, with the positive direction always toward the outer bank; and ω is the boundary inclination angle, defined according to: The calculation of q rt requires Equations 14-20: where K B is the dimensionless Brooks coefficient (Brooks & Shukry, 1963); τ c0 is the critical shear stress on a horizontal bed; μ s is the static friction coefficient; and θ rp is the particle angle of repose (Wiberg & Smith, 1987). A summary of the derivation of K B is provided in Vanoni (1975). When using the assumption of transverse slope only (longitudinal slope negligible), the relation of Kovacs and Parker (1994;Equation 18 therein) is identical to the above formulation. When δ = 0, K B is equivalent to the more commonly known expression from Lane (1955; Equation 1 therein). The bedload transport rate vector magnitude is calculated according to a general Meyer-Peter and Müller (1948) type formulation: where q* is the dimensionless volumetric bedload transport rate per unit width; g is the gravitational acceleration constant; R s is the sediment submerged specific gravity; D 50 is the median sediment grain size, which in the present analysis is the specified uniform grain size; K 0 is an empirical constant, which equals 5.7 in Fernandez Luque and Van Beek (1976); τ* is the Shields number; and τ* c is τ c made dimensionless using Equation 18. The bedload transport rate vector direction is calculated according to: where A is a coefficient that is assumed constant in linear formulations (e.g., Parker, 1984) and has dependence on β, ω, and δ in nonlinear formulations that account for steep slopes (Waterman & García, 2019). Equation 19 indicates that the direction of the bedload transport rate vector is dependent on the boundary shear stress vector direction and the downslope gravitational force. A comprehensive analysis of the past usage of Equation 19 is provided in Baar et al. (2018). Following calculation of |q*| and β, q rt is calculated as: Direct calculation of q rt is possible if values for τ b (both magnitude and direction) and ω are known, which highlights the importance of being able to specify ω at the location where r A will be chosen. In previous analyses, Hasegawa (1989) and Parker et al. (2011) specified ω as a constant value over the entire bank and set the position r A for usage in their integral formulations to the thalweg position; that is, the inner boundary of the bank region. The thalweg is the most logical location for r A , as it is a nonarbitrary feature that can be identified on any cross-section. However, the value of ω should not be assigned arbitrarily, as it clearly influences q rt , through both τ c (Equations 14 and 15) and β (Equation 19). Reinforcing this point, Hasegawa (1989) considered the tanω term in Equation 19 to be the primary driver in the decomposition of q; more specifically, he neglected the tanδ term in Equation 19 as being of a lower order of magnitude than the tanω term. This analysis reexamines that assumption by analytically establishing the value of ω at the thalweg.

Transverse Slope at the Thalweg
As the thalweg is defined as the position of minimum η on the cross-section, ∂η/∂r would equal zero if the function of η versus r was a continuous, differentiable function; thus, ω would equal zero. However, it is not evident that η versus r is necessarily a continuous, differentiable function; an abrupt grade-break at the thalweg is the alternative, which was used by Parker et al. (2011;Figure 10 therein). Note that the Parker et al. (2011) cross-section was specifically intended not to be a parallel retreat cross-section, as it allowed the region bankward of the thalweg and barward of the thalweg to adjust independently.
Two different conditions may prevail at a potential grade-break: first is that in which q rt varies continuously and smoothly with respect to r; second is that in which an abrupt discontinuity in the relationship between q rt and r exists at the thalweg, referred to as a shock in Parker et al. (2011). In the first scenario, an abrupt topographic break is not possible. Continuous variation of q rt with respect to r requires that q rt | rthal , calculated from both ω values bounding the thalweg, be equal. For clarity, these ω values are referred to as ω left and ω right . To ensure sediment is evacuated from the bank region for a migrating cross-section, this q rt | rthal value must be negative.
Although ω may have multiple values at a point in space, as in the case of the grade-break, τ* must be single-valued at a point in space on physical grounds. As helical flow causes τ* to be directed inward (to the left, negative δ) at a bend, bedload left of the thalweg is driven upslope, and bedload right of the thalweg is driven downslope; intuitively, q rt should be more negative when comparing any ω right to ω left , and this is indeed the case, based on the following considerations. Noting that q rt | rthal equals |q|sinβ, the transverse slope influences both |q| and β. Right of the thalweg, an increase in inclination tends to increase |q| through a decrease in τ* c (Equations 14 and 15); while an increase in inclination makes β more negative (Equation 19). Thus, q rt monotonically becomes more negative with increased inclination. Left of the thalweg, the situation is more complex. An increase in inclination causes β to become more positive (Equation 19), but |q| does not monotonically increase or decrease, due to the non-monotonic behavior of τ c when τ* has an upslope component. Equations 14 and 15 indicate that with increased inclination τ* c initially increases to a maximum greater than τ* c0 and then decreases. While a mathematically elegant solution was not obtained to demonstrate that q rt | rthal cannot be the same for different values of ω under the specific conditions evaluated, the following expression obtained from Equations 14-20 is used for the demonstration: The variables K B and β are each a function of ω and δ, which are calculated from Equations 15 and 19, respectively; |τ*| and δ are treated as parameters and are evaluated across a reasonable parameter space based on gravelbed streams. The variable ω right is set to 0, which yields the least negative value of q rt | rthal that can be obtained for ω right . For each combination of |τ*| and δ evaluated, ω left is evaluated over the entire range from 0° to −40° at 1° intervals to find the maximum value of the ratio; a ratio of 1 for any ω left < 0 indicates the grade-break is 9 of 33 possible. The value of |τ*| is varied from 0.046 to 0.20 in intervals of 0.002; and δ is varied from −1° to −45° in 1° intervals. All combinations of these parameters are evaluated. Constants used in the analysis are μ s = 0.84; A = 1.54; τ* c0 = 0.045. The results reveal that the maximum value of the ratio in Equation 21 is uniformly 1.0 for all combinations, associated with ω left = 0. Therefore, any ω right > 0 (resulting in a more negative q rt ) will not yield a solution; and the only solution is: In other words, the grade-break is not permissible for this first scenario that does not include a shock condition. A potential complicating factor to the above analysis is whether the parameter A used in Equation 21 is properly treated as a constant. For transverse slopes up to ∼10° relative to horizontal, this linear assumption includes minimal error compared to a more robust nonlinear formulation of A, such as that derived in Waterman and García (2019), which was based on Kovacs and Parker (1994) and Seminara et al. (2002). The nonlinear expression of A in Waterman and García (2019) is the following: where χ is a dimensionless coefficient that expresses the ratio between the mean fluid velocity impinging on a static particle versus a particle moving in the bedload layer; and μ D is the dynamic friction coefficient. Solving β in Equation 19 iteratively using the expression from Equation 24 for A, the same process of evaluating Equation 21 over the same parameter space is repeated. The results are identical; the maximum value of the ratio is 1.0 across the entire space, associated with ω left = 0.
The second scenario in which the grade-break may exist is that associated with an abrupt discontinuity in the q rt versus r relationship. Under this shock condition, the only constraints to q rt | rthal values associated ω left and ω right are that the two q rt | rthal be negative to ensure sediment evacuation from the bank region; and the integral condition (Equation 12) must have the same solution for both the region bankward of the thalweg (in which r A = r thal and r B = r out ) and inward of the thalweg (in which r A = r in and r B = r thal ). With only these constraints, the potential presence of a grade-break cannot be disproven mathematically. Nevertheless, through physical reasoning, some approximate conclusions can be drawn. The τ b distribution in the vicinity of the grade-break is effectively that which exists at the corner of a trapezoidal channel. The τ b increases bankward of such a corner for a short distance, reaching its maximum value a short distance up the bank (e.g., Khodashenas & Paquier, 1999;Thorne et al., 1998), depending on the corner angle. Without any consideration of the region interior to the thalweg, for such a situation to exist under parallel retreat, Equation 9 requires that the outer bank must be erosional across its entirety; and Equation 7 requires that ∂q rt /∂r be positive over its entirety to counteract the effect of the negative q rt /r and ensure erosion. However, in the near-corner region, q rt would be more negative moving outward (∂q rt /∂r negative) where τ b increases up a slope of constant or increasing ω. This violates the requirement of Equation 9; during bank profile evolution, the near-corner region will be depositional according to Equation 7, thus tending to smooth out the discontinuity. This always occurs in the development of straight channel experiments in uniform noncohesive materials (e.g., Ikeda, 1981), and also for the few documented cross-section sets in noncohesive experimental meandering channels (e.g., Friedkin, 1945;Schumm & Khan, 1972), even though these latter cannot be presumed to constitute equilibrium parallel retreat cross-sections. While this does not rise to the level of mathematical proof, it does suggest that the formation of a grade-break is unlikely under the assumptions of this analysis, and the most justifiable shape constraint is the one developed previously, namely that ω| rthal = 0.

Migration Rate Equation
Based on the findings in the previous section, Equation 12 can be re-expressed in a form that better highlights the relevant variables that drive migration. The integral limits are set as r A = r thal and r B = r out , with the boundary condition that q rt | rout = 0. This yields the following equation: In the following subsections, justification is provided for eliminating the integral term in the numerator of Equation 25. Next, the expression of q rt | rthal is developed, which allows the migration rate to be expressed in terms of excess boundary shear stress and other dimensionless variables. Finally, the dimensionless migration rate is expressed.

Integral Term in the Numerator
Since the shape constraint does not fully specify the channel geometry through a similarity function, various terms must be left in the formulation as parameters to be specified or neglected on order-of-magnitude considerations until future advances are achieved. The integral term in the numerator of Equation 25 results from the varying shape of the control volume in the radial direction. This term can be neglected on order-of-magnitude considerations, based on the following factors. Due to the ω| rthal = 0 shape constraint at the thalweg, tanβ = tanδ from Equation 19, and the near-bed flow is toward the inside of the bend (e.g., Rozovskii, 1961). Thus, tanβ, tanδ, and q rt are negative at r thal . This requires that q rt transition from the negative value at r thal to zero at r out . Equation 10 indicates that where ∂η/∂r = 0, then ∂q rt /∂r ≈ 0; that is, it is near a local minimum in q rt versus r.
When making the conservative assumption that the large negative q rt value that prevails at the thalweg is constant over the integral limits, the integral is still smaller than the first term in the Equation 25 numerator by a factor of ln(r out /r thal ). This logarithm will generally be less than 0.1 (r out /r thal < 1.1), given that the width of the bank region is generally much smaller than the thalweg radius.

Transverse Component of Bedload Transport Rate at the Thalweg
To further refine Equation 25, q must be obtained through hydrodynamic analysis and then decomposed to express q rt | rthal . For the case of steady, uniform bend flow, the s-momentum equation for a thin control volume encompassing the full flow depth can be expressed in the following alternative forms upon neglecting the crossstream and vertical transport of s-momentum: is based on a thin control volume with vertical edges in which c 0 is a dimensionless coefficient between 0 and 1 (e.g., Ikeda & Nishimura, 1986) that accounts for momentum extraction due to fluid shear stress on the transverse face associated with the bank; h is the local flow depth; and S is the local longitudinal slope. Equation 26b is based on an irregular-shaped control volume configured such that the s-directed fluid shear stresses on the transverse faces equal 0 (e.g., Pizzuto, 1990). In this configuration, h′ is the local hydraulic radius, which is the cross-sectional area of the irregular control volume divided by its wetted perimeter. From Equations 26a and 26b: Without a similarity function to indicate the proximity of r thal to r out , c 0 must remain as a parameter. The transverse component of τ b for steady, uniform bend flow can be calculated based on the analytical treatments of van Bendegom (1947), Rozovskii (1961), and others, which take the form: where τ b,rt is the transverse component of τ b tangent to the plane of the surface; and K 1 is an order 10 dimensionless constant taking values between 7 (Engelund, 1974) and 11 (Rozovskii, 1961). The theoretical analyses that established Equation 28 were based on integration over a vertical slice of the stream unaffected by the fluid shear stress associated with the side walls. As a result, van Bendegom (1947) and Rozovskii (1961) obtained the local flow depth (h) rather than h′ in the numerator of Equation 28. However, when integrating the centripetal acceleration, the radial pressure gradient, and the boundary shear stress following van Bendegom (1947) over an irregular control volume configured to eliminate shear stress on the transverse faces, the variable that results is h′.
This is equivalent to h when the control volume is a vertical slice unaffected by side walls. The resultant boundary shear stress is used to calculate the bedload transport rate: where γ is defined to make subsequent equations more concise. The following conditions prevail at the thalweg: ω = 0; tanβ = tanδ; and τ c = τ c0 based on Equations 14 and 15. An expression for q rt | rthal is obtained using Equation 17-20 and Equation 28: where τ* s is τ b,s made dimensionless using Equation 18. Note that if γ = 1 is specified, this is equivalent to the simpler, but less rigorous, approach of calculating q s directly from τ b,s rather than from |τ b |.

Incorporating the simplifications from the previous subsections and substituting Equation 31 into Equation 25
leads to a dimensionless migration rate equation that is in the form of an excess Shields stress formulation: where M* is the dimensionless migration rate; and H nc is the thickness of the noncohesive layer with respect to the deepest point in the channel. Although additional terms, such as τ* s and γ can be more fully expressed in terms of h and r, the form of Equation 32 is preferred, as it more directly expresses the migration rate dependence in terms of the following dimensionless groupings: (D 50 /H nc ), (h/r), and τ* s . Note that when h/r = 0 (straight channel), M = 0, as the steady state is simply the threshold channel that can contain a region of active bedload transport at the thalweg. Also note that the nondimensional variables that are obtained in Equation 32 are not the only way to express the migration rate in dimensionless form. Utilizing the width B in dimensionless groupings and normalizing the length scale of the migration rate with respect to B has a long history of use (e.g., Hickin & Nanson, 1984). That would be an ideal form for the expression, but B could not be obtained from the present analytical treatment.

Modeling Framework
A simplified numerical model is presented to answer several questions that the analytical treatment cannot: (1) Is parallel retreat possible without incorporating a time-averaged cyclical process of basal erosion and mass failure within the noncohesive layer? (2) What is the resulting channel shape (η| rthal , r thal , etc.)? The first question naturally arises due to long-held misconceptions regarding the nature of fluvial erosion in noncohesive bank materials, which is rooted in the widespread usage of the linear excess shear stress formulation (Equation 1) that strictly applies only to cohesive sediment. The second question arises due to various terms being needed in Equation 32 that cannot be established without the channel shape (e.g., η| rthal and c 0 ).
A cross-sectional numerical framework very similar to that developed by Stecca et al. (2017) is implemented herein. A 1-dimensional (1D, r-direction) fixed (Eulerian) grid is used, which is discretized into numerical elements at equal spacing ∆r. The center of each numerical element is referred to as the node; calculated variables are nodal values. The constant-width region of numerical elements comprising the active channel is numbered as nodes 1 through N. An illustration of the 1D domain is provided in Figure 4. boundaries of numerical elements (i ± 1/2); the gray shading represents the active channel region of nodes that will shift rightward in time. A suitable spatial resolution is much denser than that shown; ∆r is appropriately set to <0.1 times the mean flow depth.
The active channel intermittently shifts outward in discrete single-node increments based on physical processes. Translation of the active channel region occurs due to bedload transport processes at node N. The q rt | rout = 0 boundary condition ensures that node N is degradational. The slope angle that develops between node N − 1 and node N is projected to the top of the noncohesive layer following boundary deformation as shown in Figure 5. When the projected slope intersects the top of the noncohesive layer at, or to the right of, node N + 1, then node N + 1 becomes incorporated into the active channel at the next time step. All nodes i = 1 to N + 1 are then renumbered according to: i t+∆t = i t − 1. Conceptually, this assumes the upper cohesive layer fails vertically at node N + 1 when the shift occurs. The r coordinates of nodes 1 through N maintain the same values before and after the grid translation.
Another potential cause of grid translation is mass failure within the noncohesive layer. Oversteepening of the noncohesive soil to angles greater than the friction angle (Φ) will lead to shallow planar failures (Thorne, 1982).
Although Φ (a mass property of the soil), and θ rp (a property of grains) can potentially have different values, the two are conventionally assumed to be equivalent. This is valid in the case of an unstructured noncohesive deposit, as assumed in the present analysis; for the alternative case, see Millar and Quick (1993) and Darby et al. (2007).
Oversteepening will occur when an interface develops on the slope in which τ b = τ c , causing the lower slope (τ b > τ c ) to deform, while the upper slope (τ b < τ c ) remains static. This steepens the slope near the interface and creates a topographic discontinuity (Kovacs & Parker, 1994). However, the present model simulates bankfull flow, in which the noncohesive layer is always submerged. Even where the boundary shear stress becomes very small, the transverse slope steepens to the condition in which bedload transport occurs over the entire bank profile, because as ω approaches θ rp , τ c approaches zero. Extensive simulations revealed that slopes greater than θ rp do not develop during the evolution toward the parallel retreat profile when using the assumption that Φ = θ rp ; only a specified over-steep initial condition yields a situation in which a mass failure algorithm is warranted for the noncohesive bank material. Therefore, a planar failure algorithm is not relevant to the present analysis, although details of the one developed when considering sub-bankfull flow can be found in Waterman (2017).

Model Processes
Quasi-steady, uniform, fully developed bend flow is assumed, such that all divergence terms associated with mass and momentum fluxes in the s-direction are neglected. The parallel retreat cross-section under a specified bankfull discharge (Q bf ) is sought. Two required parameters are unknown: B and η| rout . The latter cannot be set arbitrarily high on the bank, as it potentially yields a condition in which more volume is eroded outside of the thalweg than can be deposited inside of the thalweg, thus preventing the development of steady-state migration. This issue is addressed in Appendix A, along with the rationale for establishing η| rout as a function of B in Equation A12; therefore, only one unknown parameter B remains. The procedure is iterative, with an initial estimate of B and an arbitrary initial cross-sectional geometry specified. The numerical model contains two steps at each time interval, based on the quasi-steady approximation: (a) hydrodynamic calculation to obtain boundary shear stress distribution and subsequent calculation of q rt ; and (b) boundary deformation based on an explicit finite difference representation of Equation 7. For the hydrodynamic calculation, rather than solving the mass and momentum conservation equations for the given Q bf by solving the water surface elevation, ξ(r), at each time step, an alternative procedure is used to reduce computation time. The transversely sloped water surface is maintained at the same elevations for nodes 1 through N throughout all time steps. The boundary shear stress distribution is calculated at each time step based on the momentum conservation equations. Q varies in time during the period of cross-section shape adjustment, since ξ(r) is time-invariant, whereas η(r) varies in time. The hydrodynamic and boundary deformation calculations are detailed in Appendix B.
The discrete time intervals between occurrences of active channel shifts defines ensembles of changing channel geometry between shifts in the Eulerian framework; steady state is considered to have been reached when cross-sectional geometries at specified times within ensembles have elevation differences within a specified tolerance. The ensemble comparison technique is described in Appendix C. Following achievement of steady-state migration, the calculated Q from Equation B8 is evaluated according to the criteria: |Q −Q bf | ≤ 0.001Q bf . If the criterion is not satisfied, then the procedure begins again with an increased B value when Q < Q bf or a decreased B value when Q > Q bf . The final solution is insensitive to the initial cross-section, provided that some portion of the initial cross-section is deep enough to yield τ b > τ c at the bankfull flow level. All the parameters required for the model are provided in a table in the model application of Section 5.

Study Site Description
The Mackinaw River is a gravel-bed, single-thread, meandering river that is a large tributary of the Illinois River in the state of Illinois, central United States. The basin area is 2,976 km 2 (Figures 6a and 6b). The river is ∼203 km long, with a mean longitudinal slope of 0.00056 and minimal slope variation along its length. The field study site and the nearest flow gauge station are located 12.7 and 27.3 km upstream of the Illinois River confluence, respectively. The USGS (U.S. Geological Survey) gauge station, number 05568000, has been operated continuously from 1921 to 2020. The mean daily discharge is 21.9 m 3 s −1 and the 2-yr recurrence interval annual series peak flow is 237 m 3 s −1 ; an inflection in the gauge station rating curve is present at approximately Q = 241 m 3 s −1 , which suggests overbank flow is initiated at the gauge station at approximately the 2-yr recurrence flow. The bankfull flow (Q bf ) is estimated to be 240 m 3 s −1 at the field site.
The basin topography is dominated by a sequence of roughly parallel morainal ridges aligned in a general northwest to southeast direction, with the westernmost being the Shelbyville moraine (Figure 6b). The Mackinaw basin west of the Shelbyville moraine is a broad valley containing predominantly sand and gravel glacial-fluvial deposits, formed by glacial outwash torrents (Johnstone, 2003); these deposits have been reworked by the modern river and include post-glacial alluvium. The submorainal Mackinaw River located west of the Shelbyville moraine terminus is a geomorphically distinct unit of the river (Gough, 1994) subject to high migration rates, frequent cutoffs, and a relatively large width-depth ratio, due largely to the predominantly noncohesive lower bank materials. The field site is an actively migrating bend that contains composite banks with a 4 m typical height from thalweg to top of bank. The upper cohesive layer comprises approximately half of the bank height and stands nearly vertical; the lower noncohesive layer of sand and gravel comprises the remaining half of the height and is generally concave-upward in cross-section with a much milder transverse slope. As shown in Figure 6d, the study bend has migrated extensively between 1998 and 2016; the maximum bank migration at the bend apex is 102 m, or ∼5.6 m/yr. The lower bank commonly contains slump blocks from the upper layer that develop annual vegetation if allowed to persist for sufficient duration without being eroded (Figure 6e).

Field Study
The analytical and numerical treatments represent a canonical case based on assumptions that will never be realized under real conditions. Thus, the field study was performed only to ascertain whether the canonical case reasonably approximates field conditions with respect to cross-sectional shapes and migration rates predicted. Two topographic surveys were performed using conventional total station survey methods and a local horizontal/ vertical datum. The first survey was performed on 19-20 September 2013 during low flow (1.41-1.78 m 3 s −1 ). The second topographic survey was performed on 29-30 September 2015 during moderately low flow (8.16-9.10 m 3 s −1 ). Both surveys were performed after spring flood events with Q > Q bf with no geomorphologically active events between the flood events and the respective surveys, with the intent for the surveys to represent reshaped channel geometry due to recent large flood events. The flow record at the gauge station is illustrated in Figure 7.
The topography was not obtained along monumented cross-sections; rather, data points were collected to represent full coverage of the banks, channel, and point bar. The 2015 topographic survey focused on the downstream portion of the bend subjected to the highest migration rates. A triangular-irregular-network (TIN) was created for each survey, and the topography for the cross-sections were sampled from the TINs. The cross-sectional locations are shown in Figure 8. Each cross-section had similar characteristics when comparing the 2013 and 2015 data, with the sections translated outward with discernible but minor changes in shape. In Figure 9, cross-section D illustrates representative conditions; at this section, the data transects on the banks from 2013 to 2015 were directly in alignment, and thus, interpolation errors associated with sampling from the TIN in the bank region are minimized. Between cross sections C and F, the channel migrated an average of 8.15 m with a maximum of 10.8 m; the average was obtained by calculating the area between the top of bank lines shown on Figure 8 and dividing by the average length of the top of bank lines between those sections. In addition to the topographic data, surface soil samples of the bank were collected and surveyed in 2013 between cross-sections C and D. The locations of these samples are shown on Figure 9. From the water line to the top of the noncohesive layer, the soil became finer-grained upwards in a similar manner to that observed on the point bar. The grain-size distributions are illustrated in Figure 10. Discharge record from USGS station 05568000 just upstream of the site. The field surveys were performed during the periods labeled "C" and "G." Flood events are labeled "A," "B," "D," "E," and "F." Periods with no data are due to ice. The gray inset that contains the event labeled "D" includes the record from station 05568000 and the next gauge station upstream, 05567500; the latter captured more of the short winter event and is considered to be the best available approximation of discharge conditions at station 05568000.

Numerical Model Results and Comparison With Field Data
The results of three different simulations are illustrated; all parameters used in the simulations are identical, except the D 50 and τ* c0 values. Simulation 1 has D 50 = 3.5 mm and τ* c0 = 0.036. Simulation 2 has D 50 = 7.5 mm and τ* c0 = 0.040. Simulation 3 has D 50 = 15.0 mm and τ* c0 = 0.045. Few grains were recorded in the sediment samples with D > 15 mm (maximum 4% mass for grains coarser than 15 mm in sample S1), and simulations were not extended to greater D 50 than sampled in the field. The D 50 values were selected for analysis, and the associated τ* c0 values were determined according to a Shields diagram (García, 2008). The full list of parameters is provided in Table 1. Several of Table 1 parameters are introduced in Appendix B: C f is the friction coefficient, B 0 is the reference channel width in a straight reach; H 0 is the reference channel depth in a straight reach, and S c is the longitudinal slope along the channel centerline. The values for B 0 and r c were estimated based on aerial photograph analysis; H 0 was an estimate based on field observations; S c was based on longitudinal profiles of the river (Gough, 1994); and Q bf was based on data analysis from the nearby gage station. Figure 11 and Table 2 illustrate the results of the simulations.
On Figure 11, the η values are relative to a vertical datum at the floodplain surface elevation (top of cohesive layer); c 0 is back-calculated from Equation 26a using the outputs of the hydrodynamic variables. The M values are the migration rates recorded in the simulations. Using the model outputs of sediment transport rates and geometry variables confirmed that the simulated migration rates were equivalent to the calculated value using Equation 25. Simulation 1 has the sediment parameters that best represent the field data, but the result does not accurately represent the cross-sectional shape observed in the field. By increasing D 50 and its associated τ* c0 (from a   Figure 12. The vertical is exaggerated 4 times. The marks S0 through S3 are the locations of the surface soil samples. Sample S0 was from a low midchannel bar located off the cross-section; S1, S2, and S3 were sampled on the bank. Shields diagram) for Simulations 2 and 3, the resulting shapes approach the observed cross-sectional shapes. The simulation results suggest the Mackinaw River behaves according to a considerably higher τ c0 than that indicated by the grain size distribution observed in the field, which is a key issue addressed in the Section 6.
Interpretation of the relevant physical processes from the three simulations is clearest when considering the steady-state cross-section for the D 50 = 15 mm case as an initial condition, and then evaluating a change in the simulation parameters to D 50 = 7.5 mm and τ* c0 = 0.040. The excess Shields stress across the entire section is higher, as the boundary shear stress is the same as the D 50 = 15 mm simulation, but the Shields stress increases due to the reduction in D 50 . The outer bank immediately begins to migrate rapidly outward, and the greater erosion volume from the outer bank causes the thalweg to aggrade. This aggradation results in a decrease in excess Shields stress across the section, which represents one negative feedback mechanism to reduce migration rates when the migration rates are exceedingly large. As the cross-section aggrades in the deep portion that conveys the most flow, it is no longer capable of conveying Q bf , and consequently the channel width B needed to be widened. Due to the boundary condition that the transverse component of bedload transport (q rt ) is zero at the inside node, the region interior to the thalweg tends to aggrade until τ* ≈ τ* c at the inner node, which requires less flow depth for the smaller-sized sediment. The transverse slope on the point bar becomes less steep relative to the narrower channel with greater D 50 . The same general evolution occurs when starting from the steady-state cross-section associated with D 50 = 7.5 mm and modifying D 50 to 3.5 mm.
Beyond D 50 and its associated τ* c0 , the most influential parameters in the model are S 0 , C f , r c , K 0, and θ rp . To maintain the same logic as the evaluation of the D 50 influence above, the interpretation is most easily conceived beginning from the steady-state profile for D 50 = 15 mm, maintaining D 50 = 15 mm, and modifying parameters individually. A sensitivity analysis that simply compares the effect of the parameter change on M, B, r thal , η| rthal , or τ* s |r thal does not capture the ways in which these variables influence each other during shape evolution. A qualitative description is deemed more informative. An increase in S 0 increases the boundary shear stresses, which causes rapid migration and aggradation of the thalweg, similar to the effect of reducing D 50 , but the width B is reduced due to the increased velocities that cause Q to exceed Q bf . The migration rate increases substantially. An increase in C f causes negligible change in overall channel shape, thalweg elevation, or migration rate but requires an increase in channel width to convey Q bf . An increased curvature (decrease in r c ) increases the radial component of the sediment transport rate. The channel initially produces net deposition and steepening of the point bar, leading to an elevation increase at the interior boundary. This boundary change causes the channel to lose the capacity to convey Q bf and the channel consequently is required to widen. The increased elevation of the inner boundary partly offsets the steeper transverse slope on the bar, causing only a modest amount of channel deepening at the thalweg and modest increases in B and migration rate. An increase in K 0 results in negligible change in channel shape, but an increase in the migration rate. An increase in θ rp causes the outer bank to steepen, deepens the thalweg and moves it bankward, and increases the convexity of the point bar. The increased depth of the thalweg tends to increase τ* s , while the closer proximity of the thalweg to the outer bank causes a decrease in τ* s due to the increased momentum extraction. The net effect is that τ* s is reduced. The migration rate is decreased due to the increased thickness of the noncohesive layer and decreased τ* s . Figure 12 illustrates a time series of the cross-section evolution from an arbitrary initial condition to the steady-state configuration for the simulation with  Note. The Equation # column indicates the equation number in the preceding analysis and in Appendix B in which the parameter is used.

Table 1
List of Parameters Used in the Numerical Simulations D 50 = 15 mm and τ* c0 = 0.045. As the water surface elevation is kept steady in time, the discharge changes until the steady-state configuration is achieved. The initial evolution of the trapezoidal noncohesive boundary involves excavation of bed material from the bank region into a concave-upward form with concurrent deposition of material near the inner bank. As evolution progresses, the thalweg slightly over-deepens and then recovers asymptotically to the steady-state profile. Figure 13 reveals that the distribution of q rt is negative throughout the profile, reaching a maximum negative value near the thalweg and transitioning to zero toward both boundaries based on the q rt = 0 boundary conditions. Equation 10 suggests that the q rt should be near a local minimum at the thalweg, which is confirmed in the numerical model results. The integral term containing q rt included in Table 2 is in the numerator of Equation 25. In Section 3.2.1, the integral term was eliminated from further usage in the analytical treatment based on order of magnitude considerations, and Table 2 shows that the integral term is indeed of a lower order of magnitude than the other term in the numerator, q rt | rthal . Also illustrated in Figure 13 is τ* s and the reference value of the downstream component of Shields stress (τ* ref ), defined as follows: * = ℎ The τ* s calculated in the numerical model implicitly includes the effects of transverse shear stress, so that τ* s and τ* ref (which assumes no transverse shear) can be directly compared. The only significant deviation between τ* ref and τ* s occurs near the outer bank. At the thalweg, the model reveals that τ* s is only 84.8% of τ* ref . For simulations 1 and 2, in which the thalweg was more radially distant from the outer bank, the momentum extraction was less substantial, as shown by the c 0 values in Table 2.
The numerically generated cross-section for the simulation with D 50 = 15 mm and τ* c0 = 0.045 is compared to the field surveyed sections in Figure 14. As shown, the numerically generated cross-section captures the general shape of the field-surveyed cross-sections in terms of the concave upward bank region and the modest convexity on the point bar. However, the numerically generated section is generally not as deep below the top of bank and the thalweg is several meters closer in the radial direction to the outer edge of the channel.
The numerical model simulates migration under a steady Q bf . The meander bend at the field site underwent migration and cross-sectional reshaping over a range of flows that included Q bf . Therefore, comparison of the migration rates between the field site and the simulations are approximate and indicate only whether the migration rates are on the right order of magnitude. Between the 2013 and 2015 surveys, the surveyed region in Figure 8 Figure 6d shows digitized bank lines over the 18.2-yr period of migration between 1998 and 2016. The maximum migration distance was 102 m. Using the 15-min time series, Q ≥ Q bf was experienced for 47.8 days during this time period; and Q ≥ 0.8Q bf was experienced for 81.9 days. Assuming migration only occurred when Q ≥ Q bf suggests a mean migration rate at the locus of maximum migration equal to 2.1 m day −1 and assuming migration occurred throughout the period Q ≥ 0.8Q bf suggests a mean migration rate of 1.1 m day −1 . As before, the results suggest the order of magnitude predicted by the simulations is correct.

General Considerations
The analysis presented is a conceptualization intended to shed light on bank erosion and channel migration when noncohesive bank material is a dominant factor. Such a treatment highlights key physical processes and provides a basis for comparison with observations to provide clues to those factors responsible for deviations from the canonical case. The three objectives outlined in the Introduction were resolved: (a) a rational constraint on the bank shape was established (ω| rthal = 0); (b) a migration rate formulation was developed (Equation 32); and (c) parallel retreat was shown numerically to develop as conceived in the analytical treatment as a simple outward translation of the noncohesive portion of the cross-section, providing a reliable framework to integrate the cyclic process of basal erosion of the noncohesive layer and mass failure of the upper cohesive layer. The findings reveal  important differences to past treatments; they also highlight important aspects of channel shape and misunderstandings that can result from usage of the linear excess shear stress formulation (Equation 1) when modeling bank deformation in noncohesive materials.
The analytical treatment focused on the development of a steady-state migration rate equation based on integration of the Exner equation over the bank region following the logic of Hasegawa (1989), which requires calculation of the transverse component of the bedload transport rate at the thalweg (q rt | rthal ). Since q rt | rthal is dependent on the local transverse slope at the limit of integration (ω| rthal ), the slope should not be assigned arbitrarily. The ω| rthal = 0 shape constraint determined theoretically and confirmed numerically ( Figure 11) yields a substantially different migration rate equation than that of Hasegawa (1989), due to a change in the dominant factors that dictate q rt | rthal , which drives migration. Equation 19 indicates that the decomposition of the bedload transport rate vector q at any point has a contribution associated with secondary flow (the term containing tan δ) and a contribution associated with transverse slope (the term containing tan ω). Whereas Hasegawa (1989) determined that the contribution from secondary flow was small enough to be neglected, the current analysis finds that the contribution from the secondary flow is responsible for the entirety of the transverse component of the bedload transport rate at the thalweg due to the resulting horizontal transverse slope. This finding yields a migration rate (M) equation (Equation 32) that is dependent on curvature, unlike previous forms of migration rate developed using the integral treatment. Previous bank-integrated migration rate formulations implicate the helical flow associated with curvature indirectly, as it is responsible for a deeper channel near the outer bank, which results in excess velocity or excess boundary shear stress near the outer bank. The present model more directly implicates the importance of curvature in the migration process as the primary driver of the transverse rate of bedload transport due to the helical flow. The migration rate formulation developed also generalizes to the case of a straight channel, in which r = ∞ and consequently M = 0.
The analysis supports a finding of Stecca et al. (2017), who evaluated a number of variables for suitability as independent variables in lumped noncohesive bank migration formulations. Among the variables analyzed, the downstream component of sediment transport rate at the basal node was deemed among the most promising. The current analysis indicates the physical basis of this finding, in that the migration rate is based on the transverse component of the bedload transport rate. This is directly related to the downstream component of bedload transport rate through the variable β (the deviation angle of the bedload transport rate vector from the streamwise direction), which is equal to the deviation angle δ of the boundary shear stress vector at the thalweg. Another interesting finding is that the Selwyn River cross-sections they analyzed for a hydrograph that produced bankfull flow resulted in cross-sections that maintained their general shape after undergoing a migration distance greater than 15 m, similar to that documented in the present study. The shapes were also qualitatively similar to those of the present study, being concave-upward on the bank and across the thalweg, and transitioning to convex-upward on the bar. The Selwyn River is transitional between single-threaded and braided, which represents a different geomorphic setting than the purely meandering rivers treated in the present analysis. However, the physical factors dominating the migration appear to be the same for that specific case, in which helical flow through a bend caused migration of a channel comprised of noncohesive bank materials.
A substantial body of research has addressed equilibrium channel cross-sectional geometry in straight channels with noncohesive banks and equilibrium point bar geometry when the outer bank is fixed (nonmigrating). A migrating channel differs from both scenarios. The past treatments involved the channel achieving a shape such that q rt = 0 everywhere on the section. The current treatment does not invalidate those approaches; however, it does show that such a solution is specific to the case in which M = 0. The more general solution is expressed in Equation 10, which is valid for either a migrating or non-migrating steady-state condition. For M > 0, q rt will be negative across the section except near the boundaries (where q rt = 0 as specified boundary conditions). This is illustrated in the numerical results of Figure 13 but can also be substantiated theoretically. At the thalweg, a negative q rt is ensured due to the secondary flow direction in conjunction with the ω| rthal = 0 shape constraint. Assuming that the value q rt /r is small, Equation 10 indicates the negative q rt will be approximately the local minimum (∂q rt /∂r ≈ 0 where ∂η/∂r = 0). Moving in the positive r direction (bankward of the thalweg), where ∂η/∂r is positive, Equation 10 requires a less negative q rt that will transition to q rt = 0 at the boundary. Moving in the negative r direction, where ∂η/∂r is negative, Equation 10 also requires a less negative q rt that will transition to q rt = 0 at the boundary. The analysis has not advanced to the point of being able to specify η| rthal or r thal without the use of the numerical model. However, if these variables could be predicted, it is conceivable that Equation 10 could be used to calculate steady-state migrating shapes analytically in an analogous manner to past work using the q rt = 0 solution criteria (e.g., Kikkawa et al., 1976). Such an analysis would be most practicable interior to the thalweg (well away from the outer bank), where the boundary shear stress is readily determined as a function of the local flow depth (see Figure 13).
While the exact cross-section shape still cannot be specified without detailed numerical analysis of the hydrodynamics and the transverse boundary layer associated with the outer bank, both the analytical and numerical treatments indicate that the shape will have differences relative to the characteristic cross-section of Parker et al. (2011; Figure 10 therein). The convex-upward shape of the region left of the thalweg is based on the q rt = 0 solution criterion used for predicting the transverse slope on the point bar. The convex-upward shape on the bar is valid, although the results illustrated here suggest that the convexity is poorly pronounced during steady-state migration, having only a slight camber similar to the point bar platforms described by Nanson (1980), on top of which scrolls of finer-grained suspended material were deposited. Further comparison with the characteristic cross-section of Parker et al. (2011) reveals more significant differences at the thalweg and in the bank region. A concave-upwards shape develops in the bank region, with a horizontal transverse slope at the thalweg, and with increasing transverse slope moving radially outward. This is contrasted with the grade-break at the thalweg in Parker et al. (2011) in which the convex-upward topography interior to the thalweg abruptly changes to a uniform slope in the bank region. These shape differences lead to a difference in the dominant factors driving the transverse flux of sediment from the base of bank, and thus the overall migration rate. Ultimately, a typical bank profile shape will depend on the geomorphological setting, local heterogeneities, and the relative importance of the cohesive upper layer in armoring the noncohesive layer. However, the Mackinaw River field data provides one example that the cross-sectional shape suggested by this analysis is a valid form for migrating meander bends with composite banks when the noncohesive layer has sufficient thickness to not be completely covered by slump blocks.
The concave upward shape within the noncohesive material in the bank region directly contradicts the shape suggested by use of the linear excess shear stress equation (Equation 1) to model noncohesive fluvial erosion.
Applying Equation 1 to bank materials consisting of uniform noncohesive soil properties suggests the following: (a) because the boundary shear stress generally increases with depth, the erosion rate will be greatest near the base of the bank, thus, leading to steepening of the noncohesive layer in the presence of excess boundary shear stress; (b) continued steepening leads to an unstable bank slope that exceeds the soil friction angle and large-scale mass failure involving much of the bank height; (c) after failure, the basal material must be removed by fluvial action; and (d) the cycle repeats. BSTEM model results of a predominantly sand/gravel bank on the Wabash River performed by  provides a good example of the shapes resulting from (a) and (b) when applying Equation 1 to noncohesive materials. While the cyclical sequence described accurately characterizes the process for cohesive bank erosion, it inaccurately characterizes the process for noncohesive bank materials mobilized as bedload in those cases where the armoring effect of cohesive slump blocks has limited influence on the noncohesive layer. When implementing the Exner equation with a suitable bedload transport relationship, a boundary shear stress distribution that increases with depth does not cause a bank slope to steepen with depth, but rather causes the slope to relax with depth. The development of a concave-upward noncohesive bank shape in the presence of excess boundary shear stress has been well documented experimentally in straight channels (e.g., Ikeda et al., 1981). As is borne out by the present analysis, the resulting shape does not fundamentally change simply due to the addition of planform curvature. The primary difference between the straight channel case and the bend flow case is that the helical flow serves to drive a transverse component of bedload away from the base of the outer bank. This allows the bank to migrate when subjected to excess shear stress rather than aggrading the base and relaxing its slope to a threshold channel configuration as occurs in the absence of the helical flow. One way to generalize between straight and curved channels is by conceiving of a "demand" for bed material imposed by the removal of sediment at the base of the bank in the transverse direction, following Thorne (1982) regarding unimpeded removal of basal material. When the channel is straight, the bed material demanded at the base of the bank is zero; when the channel has curvature, the demand is non-zero. This demand is satisfied through the adjustment of bank shape, which affects both the boundary shear stress distribution and the transverse component of bedload. Thus, analogous to the concept of a graded channel which delicately adjusts its longitudinal slope to convey exactly the sediment fed into the stream with the available discharge (e.g., Mackin, 1948), the steady-state migrating channel delicately adjusts its transverse slope to satisfy the demand for material at the base of the bank in a manner that conforms with Equation 10.
The findings also shed light on laboratory experiments involving meandering channels in noncohesive materials. With reference to Figure 2, two primary reasons can be identified for the difficulty in achieving a condition approximating parallel retreat in laboratory experiments of meander bend evolution: (a) vegetation establishment and consequent fine sediment deposition in over-bar Region 3 is not able to suitably constrain the channel geometry on the inner bend; and (b) η| rout is commonly set too high on the outer bank relative to the elevation at which the sediment can be transported at the inside of the bend, which creates the condition in which the cross-section will be net aggradational as it migrates, as demonstrated in Appendix A. This widening and shallowing of the channel during migration ultimately devolves into a braided channel form. Both conditions could potentially be mitigated through appropriate experimental setup and procedures undertaken during an experiment to artificially incorporate point bar vegetation. However, a third issue of bank stiffening to artificially slow the outer bank migration appears to be an additional step that is warranted to best represent prototype conditions based on the present findings. As evidence, the Mackinaw River appears to behave as if the critical shear stress is considerably larger than that indicated by its grain-size distribution. This is discussed further in the following subsection.

Limitations
The highly conceptualized scenario analyzed in this study is a gross simplification of the complexities of real rivers. It is intended as a canonical case that represents one extreme of the range of possibilities that may occur with respect to the influence of the cohesive top layer on meandering rivers with composite banks; namely, the case in which the influence of the cohesive top layer has minimal effect in protecting the underlying noncohesive material through armoring. Almost innumerable factors could be cited that were neglected which may be important in certain circumstances, but several major factors stand out as warranting elaboration: (a) the mechanics and importance of the cohesive top layer; (b) the treatment of the noncohesive sediment as uniform in size; (c) transient shape changes due to variable flows and the dynamics of bar push and bank pull; and (d) the hydrodynamic treatment.
First, the effect of the cohesive top layer on bank migration processes is significant when it comprises a large percentage of the bank height. Mass failure of the material yields slump blocks that influence τ c0 of the lower bank and mitigates migration rates (e.g., Parker et al., 2011). Caution would be required when applying the present findings without field confirmation that the river qualitatively has similar cross-sectional characteristics as that found in this analysis. The findings should not be used for practical applications on streams in which slump blocks can be expected to cover the entire basal noncohesive layer. The analysis strongly suggests that additional factors other than those included in the simplified model are playing a significant role in determining cross-sectional shape and migration rate at the Mackinaw River field site. The resistance to outer bank erosion had to be artificially increased in the numerical model by increasing D 50 (and its τ c0 ) to a value considerably larger than the field-sampled D 50 to achieve a numerical cross-section that reasonably represented the field cross-sections. Field observations revealed that the majority of the noncohesive layer was free of slump blocks and the exposed sand/gravel generally contained small percentages of fines; (see Figure 10 for a representative lower bank grain size distribution). Nevertheless, residual fines would be one potential explanation for the necessity to increase outer bank resistance through increasing D 50 . Sufficient analyses have not been performed to indicate whether the reasonable calibration result obtained by increasing D 50 (and consequently τ c0 ) was simply fortunate or is more generally applicable to other rivers. To model the full range of meandering rivers with composite banks more accurately for practical application, the mechanics of the upper cohesive layer, including both fluvial erosion and mass failure, would need to be explicitly incorporated. Numerous methods are available to accomplish this, including limit equilibrium methods (e.g., Darby et al., 2002;Simon et al., 2000), stress-strain analysis (e.g., Gong et al., 2018;Samadi et al., 2013), or development of empirical functions through field studies and lab experiments (Zhao et al., 2020).
Second, given the fining upward grain size distribution observed in both the banks and the point bar at the field site, the issue of sorting and stratigraphy (e.g., Parker & Andrews, 1985) is clearly significant in determining the shape of both bar and banks. The present model assumes a single grain size; it takes no account of different τ c0 and transport deviation angles based on grain size. In the bank region, the D 50 of sediment being actively transported would vary with elevation. This would tend to cause relaxation of the transverse slope of the bank relative to that predicted for the uniform grain size, as less excess boundary shear stress (dependent on local bank inclination) would be required to satisfy the demand for bank material imposed at the thalweg. Furthermore, the concentration of large diameter grains at the base of the bank, which is common in bends, may limit the transverse component of bed material transport across the thalweg, thus slowing the migration rate due to the reduced demand for bank material. In addition to the modified τ c0 and transport deviation angles associated with grain size, the structure of the bank material is also not incorporated into the model. The banks at the field site were not comprised of loose, unstructured deposits as assumed in the model; they were structured deposits with interlocked grains that certainly influences τ c0 and θ rp (Millar & Quick, 1993). Shifting focus to the point bar, the presence of a broadly graded grain size distribution would allow the inner portion of the point bar to build to a higher elevation than would be achieved using the uniform grain size assumption. Finer sand can be transported at lower flow depths due to its smaller τ c0 , and sediment transport is required for deposition to occur. A higher elevation of the inner bar (η| rin ) requires a higher top elevation of the noncohesive layer in the outer bank (η| rout ) to achieve parallel retreat, based on the logic described in Appendix A. A higher η| rout would reduce the long-term migration rate, which is inversely proportional to the thickness of the noncohesive layer, under the assumptions of this analysis (Equation 32).
Third, the assumption of steady-state parallel retreat under bankfull flow neglects intermediate transient states.
Variable discharge will result in bank reshaping not accounted for in a lumped bank erosion model that assumes a similarity function for bank shape. Variable discharge will result in periods in which the τ b = τ c interface occurs within the noncohesive layer, leading to oversteepening and mass failure within the noncohesive layer (e.g., Nardi et al., 2012;Zhao et al., 2020). Furthermore, changes in the bank shape will yield periods in which the top of bank migrates greater and lesser distances than that suggested by the transverse flux of sediment at the base. However, over long periods of time, these transient periods of bank reshaping will average out, and the volume of sediment transported transversely across the thalweg will approach the volume of noncohesive bank material contained in the calculated migration distance. A related issue is that meander bends undergo width changes due to the dynamics of bar push and bank pull (e.g., Eke et al., 2014;Mason and Mohrig, 2019;Zhao et al., 2021). The ability to apply the steady-state approach to practical applications that involve short time scales is limited in direct proportion to the range of width change experienced by a river bend. For example, Mason and Mohrig (2019) documented width changes on the Trinity River in Texas over two time intervals within a 6-yr period. At most cross-sections, channel widening in one time interval was followed by narrowing in the next time interval or vice versa. Interestingly, bends with larger radii of curvature experienced larger normalized width changes than tighter bends. They hypothesized a closer coupling between the movement of the outer bank and the point bar in the tighter bends due to the greater strength of secondary flow. The current analysis represents full coupling between the bar and outer bank, and its dependence on bend curvature warrants additional exploration.
Fourth, the hydrodynamic treatment used in the numerical model was highly simplified. It assumed a hydrostatic pressure distribution in the vertical, which neglects the effects of vertical accelerations that are likely important close to the bank. Transverse and vertical transport of downstream momentum were neglected, and the radial component of momentum was not explicitly solved. Such simplifications were necessary to obtain boundary shear stress distributions quickly due to the large number of time steps required to achieve the steadystate migration. Applying a more sophisticated three-dimensional nonhydrostatic hydrodynamic model with the boundary shape obtained by the simplified model would be worthwhile to ascertain whether the τ b distribution (both magnitude and direction) in the near-bank region is realistic or warrants substantial improvement. A further hydrodynamic issue is that form drag was not incorporated; if form drag was significant, the τ b associated with skin friction that drives sediment motion would be reduced. Decreasing τ b to account for form drag has the same effect on migration rate as artificially increasing τ c0 as done in Section 5.

Practical Considerations
The current approach obviously and intentionally neglects important physical factors, with implications for applying the results. Application of the modeling approach developed in Section 4 outside the framework of a lumped bank erosion model in a high-resolution cross-section-based model, such as BSTEM or CONCEPTS is beyond the scope of the present analysis, due to the necessity of incorporating cohesive soil mechanics in such models. The current findings would be most useful in a lumped bank erosion model if a similarity function for the bank shape had been obtained. The analysis did not progress to the stage of being able to specify such a function. Nevertheless, the demonstration of the concave-upward bank shape with the ω| rthal = 0 shape constraint has significant modeling implications. The ω| rthal = 0 shape constraint that is valid for steady-state migration seems to reasonably represent the field data, and this constraint allows ready calculation of the transverse component of the bedload transport rate that drives the bank migration rate in either Equation 25 or Equation 32. Figure 14 shows that the field site maintained its concave upward bank shape during migration, and thus, the ω| rthal = 0 shape constraint appears reasonable even during natural (unsteady-state) migration. The numerical model whose bank erosion treatment most closely approximates the present treatment is that of Nays2DH (Shimizu et al., 2014, Equation 129 therein). In that model, a bank erosion formulation nearly equivalent to Equation 25 is used. It contains an additional term based on aggradation or degradation of the basal node. That model or any other model that explicitly calculates q at the basal node can readily incorporate the ω| rthal = 0 shape constraint to better approximate a similarity function over the bank region.
The migration rate equation (Equation 32) was theoretically conceived as an outwardly translating bend, as opposed to the more realistic scenario of a bend that undergoes some combination of concentric outward expansion, down-valley translation, and rotation. For implementation in a 2D numerical model, the migration rate represents an instantaneous local migration normal to the channel centerline at a single station, and the variables driving the migration in Equation 25 or Equation (32) would change in time and space. Differential migration rates along the bend and temporal changes in migration direction due to channel centerline evolution would allow instantaneous outward translation at a single station to be expressed in more complex forms when applied to many stations along the channel length. To account for helical flow that is not fully developed and adjusted to the local radius of curvature, the term −K 1 c o h/r in Equation 32 that originates from Equation 28 could be replaced with a local hydrodynamic calculation of tanδ. For example, the model RVR Meander (Abad & García, 2006;Motta et al., 2012) implements a spatial convolution integral of upstream curvature to calculate the local strength of secondary flow.
Models such as RVR Meander that assume a cross-sectional form rather than calculating bed evolution explicitly using the Exner equation present additional complications with respect to specifying r thal , η| rthal , and c 0 when using the migration rate equation. With the present state of knowledge, these terms cannot be predicted a priori; a user would need to develop relationships for these variables as a function of centerline radius (or δ) using either a detailed numerical model such as that in Section 4 or obtaining field data to develop reasonable relationships. The complication is that the numerical model of Section 4 is mechanistic and requires no arbitrary parameters but is not free of requiring calibration as described previously for the Mackinaw River.
Equation A3 is based on the problem definition that M and r c are constant, along with an assumption of constant λ. The functional relationship between η and r is not known a priori; but regardless of this relationship, the r values are uniformly larger in the left-hand-side than the right-hand-side of Equation A3. This requires that η| rout must be less than η| rin to achieve the equivalence since the lower limits of the integral are equal. This is demonstrated more intuitively using the approximate form of Equations A1 and A2. Equating the erosion rates and deposition rates yields: In Equation A4, the left-hand-side is necessarily less than unity, which consequently requires that η| rout be less than η| rin . For the numerical analysis of Section 5, the most physically reasonable value of η|r out is the maximum value of η| rin at which sediment can be transported and deposited, assuming that the outer bank is a previously formed inner bend lateral accretion deposit. However, Equations A3 and A4 indicate a steady-state solution does not exist for η| rout = η| rin . If η| rout is established too high on the bank, then more volume is eroded outside the thalweg than can be deposited inside the thalweg as the steady-state geometry is approached. The result is that the channel continually aggrades until no portion of the channel cross section is competent to transport sediment (M = 0), which is not a viable equilibrium geometry if the stream has any bed material sediment load. These considerations constrain the elevation of η| rout since it must be less than η| rin , while η| rin has a physical limit based on the maximum elevation at which bankfull flow is competent to transport sediment.
The approach used herein to calculate η| rout is an approximation based on Equation A4 that is primarily intended to establish η| rout consistently over a range of conditions rather than suggesting a rigorous value of η| rout based on physical principles. The variables r thal , η| rthal , and η| rin in Equation A4 are not known a priori; in fact, the first two are key variables being sought in the numerical solution. The first of these terms is approximated as r thal = r out , and the left-hand-side of Equation A4 then simplifies to r c /(r c + B/2). The second term, η| rthal , is anticipated to be somewhat deeper relative to the floodplain surface than the reference channel depth (H 0 ) due to scour at the outside of the bend; however, the possibility exists that under certain conditions, the numerical simulation may yield a very wide, shallow channel. The simplest approximation to allow for either deepening or shallowing relative to the reference condition is to set η| rthal = −H 0 , in which the floodplain surface serves as the vertical datum. The above approximations are expressed in Equations A5 and A6.
The third unknown term in Equation A4, η| rin , can be approximated more realistically. Assuming the dominant hydrodynamic forces are associated with the downstream components of gravity and bed shear stress, the maximum possible elevation of η| rin at which sediment can be transported and deposited is approximated as follows: * Combining Equations A7-A9 yields Substituting Equations A5, A6 and A10 into Equation A4 yields: in which the floodplain surface is the vertical datum. Equation A11 yields reasonable values for η| rout with viable steady-state migration for most conditions evaluated, but values were slightly too high on the bank to prevent aggradation to a static channel for all conditions. Noting that Equation A8 may overestimate τ*| rin due to momentum extraction at the side boundary not accounted for in the equation and that a transverse slope may warrant a somewhat greater value than τ* c0 in Equation A7, a coefficient was added to the right-hand-side of Equation A8. Through trial and error, the coefficient 1.2 yielded a steady-state solution over all the evaluated conditions. The modified version of Equation A11 used to set η| rout in the simulations is thus:

Appendix B: Hydrodynamics and Boundary Deformation Treatment
The key requirement of the hydrodynamic model is to reasonably represent the boundary shear stress distribution. The essence of the present treatment is the following assumptions: (a) the s-momentum equation is dominated by the driving force due to gravity and the resisting force due to fluid shear stresses; and (b) the r-momentum equation solution can be represented in an approximate manner using simplified formulations for transverse water surface slope and the deviation angle (δ) of the bed shear stress vector with respect to the s-direction.
The velocity vector is represented by three components u, v, and w, which represent the streamwise, radial, and vertical components, respectively. In the near-bank region where the w component of velocity is non-negligible, a three-dimensional hydrodynamic model accounting for non-hydrostatic pressure would clearly be the most accurate treatment; however, it comes at a significant computational expense. Most morphodynamics models use depth-averaged hydrodynamics, with the assumption of hydrostatic pressure. In the modified cylindrical coordinates for steady, uniform bend flow with hydrostatic pressure assumption, the mass and momentum (s and r) conservation equations before depth-averaging are as follows: where g is the gravitational acceleration constant; ν T is the turbulent eddy viscosity; and ν is the molecular kinematic viscosity of water. The z momentum equation, not shown, simplifies to a linear relation between pressure and depth below the free surface.
The present model uses a simplified treatment of the momentum equations. The s-momentum treatment is described first. Neglecting terms such as the radial and vertical transport of downstream momentum represents a somewhat crude approximation of the s-momentum equation when treating the near-bank region in bend flow. However, the simplified treatment that includes only gravity and boundary shear stress has proven to yield reasonably accurate results in meandering river morphodynamics simulations focused on bank erosion Motta et al., 2014). The fluid shear stress on transverse faces of a control volume is considered a potentially dominant term in the s-momentum equation, particularly in the near-bank region. Therefore, control volumes are configured such that no s-directed shear stresses exist on the transverse faces (e.g., Pizzuto, 1990); an equally valid alternative is to establish control volumes configured vertically with transverse shear explicitly incorporated (e.g., Xu et al., 2019). The nonvertical control volumes used in the hydrodynamic treatment warrants the clarification that the nodal positions in the 1D domain represent positions on the bottom boundary.
The specific method implemented to configure the control volumes with no fluid shear stresses on the transverse faces is the merged perpendicular method of Khodashenas and Paquier (1999), which is capable of treating concave-upwards, convex-upwards, and sharp-angled boundaries in the r-z plane. In this method, rays are extended upward-normal from the channel boundary until intersecting with the water surface or an adjacent ray. When two adjacent rays intersect, they are merged into a single ray segment beyond the intersection point; the merged segment is deflected intermediate to the orientations of the original rays. This yields a subdivision of the cross-section into adjacent nonoverlapping polygonal control volumes bounded by the rays. The rays originate from the bottom boundary at the lateral edges of each numerical element (i ± ½; see Figure 4). An illustration of the irregular control volumes based on this method are provided in Figure B1. Rays are projected perpendicular to the boundary angle ω; for the ray projection, ω i+1/2 is calculated as tan −1 ((η i+3 − η i-2 )/(r i+3 − r i-2 )), where the subscripts indicate the nodal position. Calculating ω i+1/2 in this manner effectively smooth the surface and is required to prevent numerical instability. The appropriate number of nodes to include in the smoothing decreases as the user-specified ∆r value increases. Near nodes 1 and N, the number of nodes incorporated in the ω i±1/2 calculation is decreased to keep the smoothing window from extending beyond the limits of the active channel. Once the polygonal control volumes are geometrically defined, the local τ b,s is calculated based on the force balance applied to the control volume. The downstream component of the water surface elevation gradient in fully developed bend flow is: where S c is the longitudinal slope along the channel centerline. The s-dimension at the bottom boundary equals r∆θ. The uniform ∆r spacing yields wetted perimeters of the numerical control volumes equal to ∆r/cosω. For all the calculations following the ray projection, tanω i is calculated using central differencing between nodes i + 1 and i − 1, except at nodes 1 and N, which are one-sided calculations. Balancing the gravitational driving force with the boundary shear force acting on a polygonal control volume yields the following: where A i represents the cross-sectional area of the polygon in the r-z plane. An assumption of Equation B5 is that the s-dimension and radial coordinate of the control volume are approximately constant over the entire height of the polygon, which is reasonable, provided that the ray length is small relative to r. Following calculation of τ b,s for each numerical node, a boxcar smoothing procedure is implemented that averages the boundary shear stress at the numerical node with the three numerical nodes on each side of it. The smoothing procedure has the same effect as the approach recommended by Khodashenas and Paquier (1999) that involves reducing the number of numerical elements for the boundary shear stress calculation; however, the smoothing eliminates the need to subsequently reassign intermediate τ b,s values to the more densely spaced numerical elements through a curve-fitting procedure. Otherwise, the numerical routine is exactly as outlined in Khodashenas and Paquier (1999). Following calculation of τ b,s , the polygon-averaged velocity (U i ) is calculated using a standard bed shear stress closure, followed by the total discharge: where C f is the dimensionless friction coefficient, which is treated as a constant parameter in this analysis.

Treatment of the r-momentum equation, Equation B3
, is based on van Bendegom (1947) and Rozovskii (1961) as described in Jansen et al. (1979;Ch. 2.2.8 therein). Such analyses are based on fully developed bend flow in the region away from the bank, which allows usage of the hydrostatic pressure assumption and neglecting all terms in Equation B3 except the following: the first term (pressure gradient), the second term (r-shear on z-faces), and the last term (centripetal acceleration). In the near-bank region, where v and w will be of the same order of magnitude, the error introduced through this simplification must be considered significant; however, in the present state of knowledge, a tractable method to integrate the equations to yield suitably simple formulations for use in morphodynamics models does not exist. The goal is to establish reasonable approximations of the transverse water surface slope and τ b,rt . The Jansen et al. (1979) analysis based on Rozovskii (1961) is used for the transverse water surface slope: where α 1 is a dimensionless coefficient whose value is very close to 1.0; and κ (=0.41) is von Karman's constant. Although dξ/dr is not strictly constant across the section, a representative value is sought for simplicity. A characteristic value of U is the depth-averaged velocity of a reference straight reach of the channel (U 0 ), characterized by a mean depth (H 0 ), and a roughness coefficient (C f0 = C f ). The depth-averaged s-momentum equation for the reference straight reach, assumed to have much greater width than the depth, allows U 0 to be calculated as: Equation B12 provides a check on the reasonableness of the S c and C f parameters selected: Using U 0 and r c as representative values in Equation B9 yields: The current approach holds ξ| rout constant at the outer bank floodplain elevation throughout the simulation; the remainder of the water surface is calculated according to: At each time step, following calculation of the τ b,s distribution, τ b,rt is calculated using Equation 28 and |τ b | is calculated according to Equations 29 and 30.
Following the hydrodynamic calculations, the boundary deformation calculations are made. The magnitude of the bedload transport rate vector q at each node is calculated using Equations 14, 15, and 17. The direction of q is calculated using Equations 19 and 24; since Equation 19 is an implicit equation when using the nonlinear form of A in Equation 24, a bisection scheme is used to solve tanβ numerically. The q rt component of q is calculated according to Equation 20. Following calculation of q rt at all nodes, a temporally and spatially discretized form of Equation 7 is used to calculate the bed elevation change at each time step: where the subscripts after variable names indicate the spatial position and the superscripts indicate the time step. The formulations for ∆q rt /∆r at nodes 1 and N are based on the boundary conditions that 1−1∕2 = 0 and +1∕2 = 0 . Each represents a zero flux into one boundary of the numerical element, while the flux across the other boundary is approximated as an average q rt between that node and the adjacent node within the domain. The interior node formulation is based on central differencing. The time-stepping represented in Equations B15 and B16 is a forward time centered space (FTCS) numerical scheme, which is conditionally stable for diffusive processes, but is unstable for advective processes unless smoothing is incorporated. The smoothing implemented to develop the boundary shear stress distribution and achieve stability is described previously. An upwinding scheme applied to Equations B15 and B16 was attempted to overcome the need for smoothing but resulted in more pervasive instability issues.
where the subscripts indicate the following: i, the spatial node; j, the cross-section evaluated within the ensemble; Ens1, Ensemble 1; and Ens2, Ensemble 2. The error term E 1 is the root mean square deviation over all nodes and cross-sections evaluated, and E 2 is the maximum individual deviation at a node. The steady-state criteria used requires both of the following conditions be satisfied: (a) E 1 ≤ 1 × 10 −5 H 0 ; and (b) E 2 ≤ 1 × 10 −4 H 0 .  Volumetric bedload transport rate per unit width vector [L 2 T −1 ] q* Dimensionless volumetric bedload transport rate per unit width [-]