Use of Multidimensional Models to Investigate Boundary Shear Stress through Meandering River Channels

Three-dimensional hydraulics were simulated through a wide range of synthetically generated meandering river channels to determine how channel curvature and width would correlate with the maximum boundary shear stress. Multidimensional models were applied, similar to a computational flume to simulate a wide range of 72 meandering channels, developed from sine-generated curves. Cannel sinuosity ranged from 1.1 to 3.0 and included five consecutive meander bends. Longitudinal slopes of the various channels spanned four orders of magnitude, while bankfull discharges spanned three orders of magnitude. Using results from one-half of the simulation sets, an empirical correlation was found to predict the maximum boundary shear stress as a function of dimensionless ratios of channel curvature and width. The remaining simulation sets were used for verification. Multidimensional models were used to simulate channel hydraulics to efficiently investigate a wide range of channel sinuosity, width/depth ratios, bankfull discharges, and valley slopes. When simulating such a wide range of channel conditions, multidimensional models offer a more efficiency method of generating consistent datasets than either field studies or physical modeling. This paper demonstrates how multidimensional models can be used to identify important hydraulic relationships that are otherwise difficult to determine.


Introduction
Meandering river channels naturally migrate across their floodplains and sometimes erode terrace banks at the edge of floodplain boundaries [1][2][3][4][5][6]. This migration process can become especially important to people living near a river and to agencies planning or maintaining infrastructure along a river [7]. An improved understanding of river hydraulics and the resulting boundary shear stress that can initiate channel bank erosion and meander bend migration is needed.
The objective of the research described in this paper is to better understand the key variables that predict increased boundary shear stress caused by curvature of meandering river channels of the type found in nature. Multidimensional models were applied to a wide range of synthetic meandering river channels to develop carefully controlled datasets. Then, the datasets were evaluated to develop empirical relationships that estimate the increased boundary shear stress.
Eventually, the empirical relationships could be used to aid in streambank protection design or help estimate the potential rate of meander bend migration. Rather than trying to obtain measured data from a wide range of natural channels, or from a wide range of laboratory flume experiments, the unique aspect of this research was the application of multidimensional models to develop carefully controlled datasets representing meandering river channels.
A meandering river channel is characterized by a sequence of smooth bends ( Figure 1). Erosion tends to occur along the outer bank, forming pools, with deposition of point bars along the inner bank. Erosion along the outer bank and deposition along the inside bank leads to the migration of meander bends both across the floodplain and along the downstream valley axis. Meandering river channels have been observed throughout the world at all scales, all continents, and latitudes, in ice [1,4,5], and even on the planet Mars [8] and Saturn's moon Titan [9]. The meander bend migration process has been observed in the field by Leopold et al. [1], Hickin [10], Hickin and Nanson [11], Hickin and Nanson [12], Parker [13], as well as in the laboratory by Friedkin [14] and Chang et al. [15]. Natural geomorphic features of a river channel, including meander bends, result from the continuous and dynamic interaction between the sediment-carrying streamflow and the erodible banks of the channel [16]. The major factors affecting alluvial stream channel forms have been listed by Lagasse et al. [17] as follows: • Stream discharge (magnitude, duration, and frequency), temperature, and viscosity; • Sediment load and grain size; The width, depth, and alignment of a meandering river channel continually adjust over time in response to changing upstream water discharge and sediment load and within the constraints of the valley slope, local geology, and human structures [18][19][20]. Yang [2] attributes river meandering to the minimum rate of energy dissipation, subject to the constraints applied to the system. Natural geomorphic features of a river channel, including meander bends, result from the continuous and dynamic interaction between the sediment-carrying streamflow and the erodible banks of the channel [16]. The major factors affecting alluvial stream channel forms have been listed by Lagasse et al. [17] as follows: • Stream discharge (magnitude, duration, and frequency), temperature, and viscosity; • Sediment load and grain size; The width, depth, and alignment of a meandering river channel continually adjust over time in response to changing upstream water discharge and sediment load and within the constraints of the Water 2020, 12, 3506 3 of 76 valley slope, local geology, and human structures [18][19][20]. Yang [2] attributes river meandering to the minimum rate of energy dissipation, subject to the constraints applied to the system.
However, the precise cause of channel meandering remains undefined. Ashworth, et al. [21] suggested that the process began with coherent flow structures that were in the order of a channel width. Large eddies produced in the flow field of a straight channel induce a sinuous path in the streamlines of maximum velocity. Seminara [16] used a linear model of flow and bed topography and showed that meanders behave as linear oscillators, which may resonate at some distinct values of the channel width-to-depth ratio and of the ratio of channel width to meander wavelength. Resonance excites a natural mode of oscillation in the form of stationary alternate bars. When the resonance barrier is crossed, the direction of meander migration is reversed.
As streamflow is conveyed through a meander bend, the centrifugal force causes the water surface elevation to be slightly higher on the outside of the bend (Figure 2). This super elevation of the water surface causes a pressure imbalance perpendicular to the main flow direction, which then causes a secondary flow current [16,22]. The secondary current is in a downstream helical motion, which is downward near the outer bank, back across the channel bottom toward the shallow inner bank, and then back across the channel near the water surface.
Water 2020, 12, x FOR PEER REVIEW 3 of 76 However, the precise cause of channel meandering remains undefined. Ashworth, et al. [21] suggested that the process began with coherent flow structures that were in the order of a channel width. Large eddies produced in the flow field of a straight channel induce a sinuous path in the streamlines of maximum velocity. Seminara [16] used a linear model of flow and bed topography and showed that meanders behave as linear oscillators, which may resonate at some distinct values of the channel width-to-depth ratio and of the ratio of channel width to meander wavelength. Resonance excites a natural mode of oscillation in the form of stationary alternate bars. When the resonance barrier is crossed, the direction of meander migration is reversed.
As streamflow is conveyed through a meander bend, the centrifugal force causes the water surface elevation to be slightly higher on the outside of the bend (Figure 2). This super elevation of the water surface causes a pressure imbalance perpendicular to the main flow direction, which then causes a secondary flow current [16,22]. The secondary current is in a downstream helical motion, which is downward near the outer bank, back across the channel bottom toward the shallow inner bank, and then back across the channel near the water surface. In river channels with non-erodible banks, flow along the inner bend accelerates relative to the outer bend. Farther downstream, the secondary flow transfers momentum towards the outer bend, moving the thread of high velocity from the inner bank to the outer bank [16]. The establishment of a bar-pool pattern creates an additional component of the secondary flow, which further modifies the bed topography.
If the channel bed is erodible, asymmetry in the velocity and boundary shear stress distributions rapidly leads to the generation of pools, riffles, and alternate bars [23]. There is a lag effect of the hydraulics through a meander bend where the curvature effects develop some distance downstream. This distance is referred to as the planform phase lag and can lead to downstream migration of the meander bed ( Figure 3). z L R c Figure 2. Typical river channel cross section with super elevation around a meander bend that causes a pressure imbalance and induces a secondary flow current. The arrows denote the direction of downstream flow on the water surface (light blue shading) and secondary circulation in cross section (medium blue shading).
In river channels with non-erodible banks, flow along the inner bend accelerates relative to the outer bend. Farther downstream, the secondary flow transfers momentum towards the outer bend, moving the thread of high velocity from the inner bank to the outer bank [16]. The establishment of a bar-pool pattern creates an additional component of the secondary flow, which further modifies the bed topography.
If the channel bed is erodible, asymmetry in the velocity and boundary shear stress distributions rapidly leads to the generation of pools, riffles, and alternate bars [23]. There is a lag effect of the hydraulics through a meander bend where the curvature effects develop some distance downstream. This distance is referred to as the planform phase lag and can lead to downstream migration of the meander bed ( Figure 3).  Cohesive and erodible channel banks are a necessary condition for the initiation of a meander bend, although bank erosion is not the direct cause of meandering [24]. The presence of deep-rooted vegetation is important for providing bank cohesion which constrains channel migration and allows deeper and narrower channels to develop [16]. The sinuous path of the maximum velocity streamlines initiates a pattern of bank erosion that marks the onset of channel meandering [14,25,26]. Where the streamline with maximum velocity converges with the channel bank, local scour of the channel bed undercuts and erodes the bank and promotes bank retreat and meander bend migration [27,28]. The rate of sediment deposition along the inside channel bend depends on the upstream sediment supply. Channel width can remain relatively constant over time only if the rates of bank erosion and point bar deposition are in balance. Brice [29] discovered that channels with constant width were relatively stable, whereas channels that were wider at bends were more dynamic. High sinuosity and equal-width streams were the most stable, whereas other equal-width streams of lower sinuosity were less stable, and wide bend streams had the highest erosion rates.
Seminara [16] argued that meander bend formation must be investigated through erosional mechanisms driven by the three-dimensional structure of the flow field. River channel morphology and rate of evolution are strongly determined by the rates that water discharge and sediment load are supplied from upstream [17]. The local distribution of velocity and shear stress and the bed and bank sediment characteristics control channel behavior. The direction and rate of channel migration are a function of the following channel characteristics [17]: local curvature and sinuosity, channel shape (width/depth ratio), longitudinal slope, discharge, and bank strength.
Hooke [30] described the distribution of sediment transport and shear stresses in a meander bend. The bank erosion equation by Ikeda et al. [31] led to empirical and theoretical studies to predict meander bend migration. Parker [32] and Odgaard [33,34] developed meander migration models by utilizing similar approaches as Ikeda et al. [31]. Smith and Mclean [35] and Mclean and Smith [36] were the first to develop dynamically correct two-dimensional models and explain the mechanics of flow through a meander bend. Nelson and Smith [37] combined a bed-load transport algorithm with the numerical model for flow and boundary shear stress fields to investigate point bars in meandering channels and alternate bars in straight channels.
Johannesson and Parker [38] developed an analytical model to simulate the increase in channel velocity through a meander bend based on a linearization of the two-dimensional equations for depth-averaged flow and a sediment transport equation to predict the cross-channel slope. The flow variables were the sum of the following two parts: (1) flow in straight channel and (2) flow deviation due to a slightly curved channel. Then, the flow variables were substituted into three-dimensional (3D) flow equations. The equations were then simplified and grouped into the terms responsible for the straight channel solution and those due to the channel curvature. The equations became ordinary Cohesive and erodible channel banks are a necessary condition for the initiation of a meander bend, although bank erosion is not the direct cause of meandering [24]. The presence of deep-rooted vegetation is important for providing bank cohesion which constrains channel migration and allows deeper and narrower channels to develop [16]. The sinuous path of the maximum velocity streamlines initiates a pattern of bank erosion that marks the onset of channel meandering [14,25,26]. Where the streamline with maximum velocity converges with the channel bank, local scour of the channel bed undercuts and erodes the bank and promotes bank retreat and meander bend migration [27,28]. The rate of sediment deposition along the inside channel bend depends on the upstream sediment supply. Channel width can remain relatively constant over time only if the rates of bank erosion and point bar deposition are in balance. Brice [29] discovered that channels with constant width were relatively stable, whereas channels that were wider at bends were more dynamic. High sinuosity and equal-width streams were the most stable, whereas other equal-width streams of lower sinuosity were less stable, and wide bend streams had the highest erosion rates.
Seminara [16] argued that meander bend formation must be investigated through erosional mechanisms driven by the three-dimensional structure of the flow field. River channel morphology and rate of evolution are strongly determined by the rates that water discharge and sediment load are supplied from upstream [17]. The local distribution of velocity and shear stress and the bed and bank sediment characteristics control channel behavior. The direction and rate of channel migration are a function of the following channel characteristics [17]: local curvature and sinuosity, channel shape (width/depth ratio), longitudinal slope, discharge, and bank strength.
Hooke [30] described the distribution of sediment transport and shear stresses in a meander bend. The bank erosion equation by Ikeda et al. [31] led to empirical and theoretical studies to predict meander bend migration. Parker [32] and Odgaard [33,34] developed meander migration models by utilizing similar approaches as Ikeda et al. [31]. Smith and Mclean [35] and Mclean and Smith [36] were the first to develop dynamically correct two-dimensional models and explain the mechanics of flow through a meander bend. Nelson and Smith [37] combined a bed-load transport algorithm with the numerical model for flow and boundary shear stress fields to investigate point bars in meandering channels and alternate bars in straight channels.
Johannesson and Parker [38] developed an analytical model to simulate the increase in channel velocity through a meander bend based on a linearization of the two-dimensional equations for depth-averaged flow and a sediment transport equation to predict the cross-channel slope. The flow variables were the sum of the following two parts: (1) flow in straight channel and (2) flow deviation due to a slightly curved channel. Then, the flow variables were substituted into three-dimensional (3D) flow equations. The equations were then simplified and grouped into the terms responsible for the straight channel solution and those due to the channel curvature. The equations became ordinary differential equations that could be solved analytically or through relatively simple numerical methods. Sediment transport is assumed to be a function of the local velocity and shear stress and is used to compute the cross-channel bottom slope. Sun et al. [39,40], rederived the equations to account for multiple sediment grain sizes.
Cherry [41] concluded that the accuracy and applicability of channel migration models were limited by the assumptions of constant discharge and channel width. Cherry [41] also noted that the local factors that influence the rate of bank erosion required further investigation and that further refinements in bend-flow modeling would not improve our predictive capability until we found a more rational way to wed the flow model to a bank erosion model.
Lagasse et al. [17] reviewed the methodology for predicting channel migration and concluded that it was difficult to predict the magnitude, direction, and rate at which channel migration would occur and that no good practical models of channel migration presently existed.
Parker et al. [42] developed a meander bend migration model where the rates of channel bank erosion and deposition were computed separately so that the channel width could widen or narrow. The eroding bank was assumed to consist of a cohesive sediment layer overlying a non-cohesive layer of sand and gravel. Bank erosion of the non-cohesive layer resulted in slump block failures of the overlaying cohesive layer. The slump blocks along the bank could slow the rate of subsequent bank erosion until the slump blocks were eroded and the bank erosion cycle continued.
Additional research has improved the understanding of bed load transport, calibration of channel roughness for numerical modeling, and the determination of environmental flows. Kuriqi et al. [43] demonstrated the effective use of physical model experiments, using a laboratory flume, to compute the critical Shields number for different hydraulic conditions and bed slope. The dimensionless Meyer-Peter and Müller number was calibrated so the computed sediment discharge matched the measured discharge. Ardıçlıoglu and Kuriqi [44] developed empirical relationships that predicted Manning's n roughness as a function of Froude number. They also found an empirical relationship between the n-value computed from measurements and the n-value needed for numerical modeling. Suwal et al. [45] evaluated the application of environmental flow assessment methods on the Kaligandaki River in Nepal and concluded that minimum stream flows of at least 30% of mean-daily discharge or those computed using annual distribution methods were more suitable for maintaining health of the riverine ecosystem.
For this research, existing two-and three-dimensional numerical models were used to simulate streamflow depth, velocity, and shear stress through a wide range of synthetic meandering river channels of the type found in nature. Sine-generated curves were used to define the meandering channel alignments with three or five consecutive meander bends. Fixed-bed, trapezoidal channels were simulated at the bankfull discharge. The model simulation results were then used to develop general relationships between dimensionless boundary shear stress and dimensionless geometric channel properties (i.e., dimensionless channel curvature and width).

Numerical Model Methods
Twelve sets of synthetic meandering channels were created for the numerical simulation of depth, velocity, and shear stress over a wide range of hydraulic conditions. Each simulation set included multiple channels representing a range of sinuosity. For each channel, the two-dimensional (2D) model, SRH-2D [46,47], was used to simulate the open channel water surface. Next, the three-dimensional (3D) model, U 2 RANS [48,49], was used to simulate streamflow velocities and boundary shear stress, using the previously simulated water surface from the 2D model.
Even though each numerical simulation took about 6 h of computer time (often overnight on a high-end laptop computer), model simulations were performed at much less time and cost than the construction and use of physical model flumes or data measurements from natural channels. The numerical model results are repeatable and there are no scaling issues. Prior to these numerous simulations, the numerical model was validated by comparison with measurements from three different physical models.

Synthetic Meandering Channels
A range of simulation sets were developed from which to develop synthetic meandering channels. The twelve model simulation sets were created based on a graph presented by Chang [50] relating bankfull discharge; channel slope, width, and depth; and median-sediment grain size for natural channels ( Figure 4 and Table 1). The graph by Chang [50] included four regions of different channel types. Model simulation sets were developed only from Region 1 (equal-width point-bar streams) and Region 3 (wide-bend point-bar streams). Region 2 represented straight braided streams, while Region 4 represented steep braided streams. Each simulation set was defined by a bankfull discharge, valley slope, and sediment grain size. Each simulation set included nine meandering channels with sinuosity of 1.10, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, and 3.00. A total of 108 synthetic channels were generated, and 72 of these channels were used for numerical model simulation.  The first three simulation sets (1, 2, and 3) contained channels with the lowest discharge and steepest slopes, while the last three sets (10, 11, and 12) included channels with the highest discharge and most mild slopes. Model results from the middle simulation sets (4 through 9) were used to develop empirical equations relating maximum near-bank shear stress with channel curvature and width. Model results from the other simulation sets (1, 2, 3, 10, 11, and 12) were used to validate the equations.
Two-thirds of the simulation sets were designed to be within Region 1 and one-third within Region 3. The eight simulation sets within Region 1 included 72 meandering channels, four of which were at the critical slope for incipient motion of sediment particles and had a sinuosity of 3.00 ( Figure 4). There were also four meandering channels at the upper boundary of Region 1 with a sinuosity of 1.1. The four simulation sets in Region 3 included 36 meandering channels, four of which were at the lower boundary of the region with a sinuosity of 3.00 and another four at the upper boundary with a sinuosity of 1.1.

Channel Dimensions
Longitudinal slopes were computed for each channel by dividing the valley slope by the sinuosity. Then, equations by Chang [50] were used to determine the surface width and depth for each synthetic channel. Region 1 characterizes channels with gentle longitudinal slopes, low velocity, small bed-material loads, and flow resistance in the lower regime of ripples and dunes [50]. Channels in this region are typically meandering and have width/depth ratios generally in the range of 4 to 20. Chang [50] provided equations to predict the channel surface width and depth in Region 1 as follows: In the channel width equation, T is the top water-surface width (m), d is the median particle grain size (mm), S is longitudinal channel slope, S c is critical slope for incipient motion of the median sediment particle size, and Q is bankfull river discharge (m 3 /s). In the depth equation, D is the center channel depth (m) at the bankfull discharge. Region 3 may include some braided river channels, but also includes sinuous point-bar rivers with riffles and pools. Chang [50] also provided equations to predict the channel surface width and depth in Region 3: The median bed-material grain size had to be assumed for each simulation set to determine the plotting position in Figure 4. According to van den Berg [51], median bed-material grain generally tends to increase with channel slope. In addition, the median bed-material grain size is most frequently between 0.1 and 1.0 mm (sand-bed channels) or between 10 and 100 mm (gravel and cobble-bed rivers). The assumed median bed-material grain sizes, for the twelve simulation sets, are presented in Table 1. The same grain size was assumed for each meandering channel of a given simulation set. Relative roughness (k s ) for the synthetic channels was computed assuming 3.5 median grain diameters (k s = 3.5 d/D) to approximate the rough boundaries of a natural channel (values ranged from 1.266 × 10 −5 to 4.988 × 10 −2 ).
The dimensions and hydraulic properties of each synthetic meandering channel are presented in Appendix A. All synthetic channels were assumed to have a constant width and a trapezoidal cross section with 2H:1V side slopes. Manning's equation [52] was used to compute the normal, bankfull flow depth and velocity as follows: In the above equations, Q is the channel discharge, A is cross-sectional area of flow, R is hydraulic radius, D is the flow depth, z is the ratio of horizontal side slope distance to unit vertical rise, B is the channel bottom width, and P is the wetted perimeter.
All flows for these simulations were subcritical with Froude numbers ranging from 0.057 to 0.539. The Reynolds numbers for these simulations ranged from 1.19 × 10 5 to 1.28 × 10 7 and in the turbulent range (dominated by inertial forces).

Meandering Channel Alignments
Langbein and Leopold [53] found that the sine-generated curve best approximates the alignment of a natural meandering river channel. Therefore, the alignments of the synthetic meandering channels were developed from sine-generated curves as: In the above equation, φ is the angle of channel alignment relative to the downstream valley axis, β is maximum angle of channel alignment, L is distance along the channel centerline, λ is meander wavelength, Ω channel sinuosity, and π is the mathematical constant (3.14286). Equation (9) produces a variable radius of curvature that is a minimum at the bend apex and increases toward infinity (at the crossing) where the meander bend ends and the curvature reverses for the start of the next downstream bend ( Figure 5).  Equation (9) was integrated to produce the following equations to compute the horizontal coordinates of the channel alignment, where x and y are the horizontal coordinates parallel and perpendicular to the downstream valley axis and l is the distance along the channel center line at positions n and n + 1: A polynomial regression equation [54] was used to relate the maximum angle of the meander bend (β) to the channel sinuosity (Ω) as follows: Ω = 0.9218 − 2.8025 + 3.3359 − 1.2941 + 1.159 (12) Each simulation set included nine synthetic stream channels with a range of sinuosity from 1.10 to 3.00. As an example, a series of meandering channel alignments were computed using a sinegenerated curve over 1.5 meander wavelengths (three meander bends) to compare the range of channel sinuosity ( Figure 6). However, the synthetic channels used in model simulations included five consecutive meander bends (2.5 meander wavelengths). Input data for the numerical models were always specified with at least three digits of precision. Equation (9) was integrated to produce the following equations to compute the horizontal coordinates of the channel alignment, where x and y are the horizontal coordinates parallel and perpendicular to the downstream valley axis and l is the distance along the channel center line at positions n and n + 1: x n+1 = x n + (l n+1 − l n ) cos(φ) (10) A polynomial regression equation [54] was used to relate the maximum angle of the meander bend (β) to the channel sinuosity (Ω) as follows: Each simulation set included nine synthetic stream channels with a range of sinuosity from 1.10 to 3.00. As an example, a series of meandering channel alignments were computed using a sine-generated curve over 1.5 meander wavelengths (three meander bends) to compare the range of channel sinuosity ( Figure 6). However, the synthetic channels used in model simulations included five consecutive meander bends (2.5 meander wavelengths). Input data for the numerical models were always specified with at least three digits of precision.
In order to keep the numerical model boundaries well away from the middle meander bends, all the synthetic channel alignments began and ended with straight channel segments (upstream and downstream from the meandering reach) ( Figure 6). Each of these straight channel segments were 10 or 15 channel widths long. This length of straight channel segments was considered long enough that the specified boundary conditions would not influence the hydraulic conditions in the middle meander bends. The hydraulic model simulation through the straight upstream channel segment was used to represent the reference bottom shear stress for the meandering channel. In order to keep the numerical model boundaries well away from the middle meander bends, all the synthetic channel alignments began and ended with straight channel segments (upstream and downstream from the meandering reach) ( Figure 6). Each of these straight channel segments were 10 or 15 channel widths long. This length of straight channel segments was considered long enough that the specified boundary conditions would not influence the hydraulic conditions in the middle meander bends. The hydraulic model simulation through the straight upstream channel segment was used to represent the reference bottom shear stress for the meandering channel.
The straight channel segments each connect to the meandering reach at the crossing points so there is a gradual transition in curvature through the bend. Many previous physical models connected a straight reach to a curve with constant radius, creating an abrupt transition.

SRH-2D Model
The SRH-2D model, developed by Lai [47], solves the two-dimensional depth-averaged Saint-Venant equations, using a finite-volume discretization. Mass conservation is satisfied both locally and globally. The model automatically simulates the wetting and drying of arbitrarily shaped mesh cells. The model can simulate both steady and unsteady flows as well as subcritical, trans-critical, and supercritical flow. The Manning's roughness coefficient is used to represent flow resistance and is usually the only calibration parameter. Lai [47]  Channel with 90° curves.
The model was validated with simulations of the Sandy River Delta Dam located near the confluence of the Sandy and Columbia Rivers, east of Portland, Oregon. Lai [47] found that quadrilateral mesh cells produced more accurate results than the same number of triangular mesh cells. Simulated water surface was less sensitive to the choice of mesh types or size, while the velocity distribution was more sensitive to the type and number of mesh cells. The parabolic turbulence model was found to be adequate for the simulation of water surface elevation and velocity distribution for The straight channel segments each connect to the meandering reach at the crossing points so there is a gradual transition in curvature through the bend. Many previous physical models connected a straight reach to a curve with constant radius, creating an abrupt transition.

SRH-2D Model
The SRH-2D model, developed by Lai [47], solves the two-dimensional depth-averaged Saint-Venant equations, using a finite-volume discretization. Mass conservation is satisfied both locally and globally. The model automatically simulates the wetting and drying of arbitrarily shaped mesh cells. The model can simulate both steady and unsteady flows as well as subcritical, trans-critical, and supercritical flow. The Manning's roughness coefficient is used to represent flow resistance and is usually the only calibration parameter. Lai [47] verified the SRH-2D model by simulating the following types of open channel flows: • Subcritical flow in a 1D channel; • Transcritical flow in a 1D channel; Channel with 90 • curves.
The model was validated with simulations of the Sandy River Delta Dam located near the confluence of the Sandy and Columbia Rivers, east of Portland, Oregon. Lai [47] found that quadrilateral mesh cells produced more accurate results than the same number of triangular mesh cells. Simulated water surface was less sensitive to the choice of mesh types or size, while the velocity distribution was more sensitive to the type and number of mesh cells. The parabolic turbulence model was found to be adequate for the simulation of water surface elevation and velocity distribution for most flow areas without flow separation. However, the k−ε turbulence model was recommended for cases of flow separation or steep banks. Use of measured water surface elevations was recommended for calibrating the Manning's roughness coefficient. The SRH-2D model has been widely used and endorsed by the U.S. Federal Highway Administration for 2D hydraulic analysis [55].

U 2 RANS Model
The 3D model, U 2 RANS, was developed by Lai et al. [48,49] as an implicit, finite volume model. The authors applied and validated the model for a wide variety of practical problems. The model solved the three-dimensional Reynolds-averaged, Navier-Stokes, turbulent flow equations to simulate the distributions of flow velocity, static pressure, and shear stress for each mesh cell. Like the SRH-2D Model, mass conservation in the U 2 RANS model was satisfied both locally and globally and the numerical formulation was applicable to arbitrarily shaped cells. The model solved the Reynolds-averaged Navier-Stokes equations for steady, incompressible, turbulent flows using the standard k-ε turbulence model. The governing equations were developed for non-hydrostatic pressure distribution.
Lai [48] verified and validated velocity distributions using the U 2 RANS model by simulating turbulent flow through an S-shaped meandering channel with a trapezoidal cross section. This S-shaped channel consisted of two identical 90 • curves (right turn then left turn), connected by a short straight channel segment. Simulation results compared very well with experimental data from a physical model flume. Lai [49] additionally validated the model velocity distributions by simulating flows in a hydro turbine draft tube and in the forebay of Rocky Reach Dam, WA, USA. Simulation results were found to agree very well with laboratory and field measurements.

Additional Model Velocity Validation
For additional validation of the model velocity simulations, the SRH-2D and U 2 RANS models were applied to simulate a large-scale physical model (laboratory flume) at Colorado State University (CSU) with two channel curves [7].
The physical model channel was 57.5 m long and included two curves of different, but constant radii and width, separated by a short transitional straight reach [56,57] (Figure 7). The upstream channel curve was to the right, while the downstream curve was to the left. The channel width of the downstream curve was about three-fourths the width of the upstream curve. and the numerical formulation was applicable to arbitrarily shaped cells. The model solved the Reynolds-averaged Navier-Stokes equations for steady, incompressible, turbulent flows using the standard k-ε turbulence model. The governing equations were developed for non-hydrostatic pressure distribution.
Lai [48] verified and validated velocity distributions using the U 2 RANS model by simulating turbulent flow through an S-shaped meandering channel with a trapezoidal cross section. This Sshaped channel consisted of two identical 90° curves (right turn then left turn), connected by a short straight channel segment. Simulation results compared very well with experimental data from a physical model flume. Lai [49] additionally validated the model velocity distributions by simulating flows in a hydro turbine draft tube and in the forebay of Rocky Reach Dam, WA, USA. Simulation results were found to agree very well with laboratory and field measurements.

Additional Model Velocity Validation
For additional validation of the model velocity simulations, the SRH-2D and U 2 RANS models were applied to simulate a large-scale physical model (laboratory flume) at Colorado State University (CSU) with two channel curves [7].
The physical model channel was 57.5 m long and included two curves of different, but constant radii and width, separated by a short transitional straight reach [56,57] (Figure 7). The upstream channel curve was to the right, while the downstream curve was to the left. The channel width of the downstream curve was about three-fourths the width of the upstream curve.
The SRH-2D model was used to simulate the channel water surface for a discharge of 0.566 m 3 /s. The simulated water surface elevations from the SRH-2D model (including super elevation along curves) were then used to define the rigid-lid elevations of 3D model grid. The U 2 RANS model was used to simulate the velocity distribution and results were compared with velocity measurements from the physical model. The upstream and downstream boundaries of the numerical models were extended ten channel widths past the boundaries of the physical model. The 3D velocity distributions across the upstream The SRH-2D model was used to simulate the channel water surface for a discharge of 0.566 m 3 /s. The simulated water surface elevations from the SRH-2D model (including super elevation along curves) were then used to define the rigid-lid elevations of 3D model grid. The U 2 RANS model was used to simulate the velocity distribution and results were compared with velocity measurements from the physical model.
The upstream and downstream boundaries of the numerical models were extended ten channel widths past the boundaries of the physical model. The 3D velocity distributions across the upstream boundary of the physical model ended up being significantly different than the simulated velocity distribution across this same location. Therefore, the simulated velocity distributions did not begin to match measured velocities until the downstream end of the first channel curve. Simulated and measured velocity distributions matched reasonably well through the second curve. This example illustrates the need for the upstream boundary of a physical model to be well upstream of the region where measurements are of great interest. The details of the 2D and 3D numerical model simulations and the comparison of simulated and measured velocity distributions are presented in Appendix B.

Model Shear Stress Validation
According to verifications by the model developer (Lai 2003a and2003b) and this research, the twoand three-dimensional numerical models (SHR-2D and U2RANS) were considered to be suitable for predicting secondary circulation and the distribution of stream velocity through channel bends. Shear stress computed by the U 2 RANS model was validated for this research by simulating flow through physical models at the University of Iowa [58] and Massachusetts Institute of Technology [59]. The 3D numerical model simulations of shear stress compared well with physical model measurements.

University of Iowa Physical Model
Lai et al. [48] previously verified the U 2 RANS model by comparing simulated water depth and velocity profiles with measurements from a physical model of open channel flow at the University of Iowa [58]. These numerical model simulations were repeated for this research to also verify the model's ability to simulate boundary shear stress. The channel of the physical model was 35.348 m long and included two identical 90 • curves of opposite direction (right curve then left curve) connected by a short straight channel segment ( Figure 8). Short, straight channel segments were also included at the upstream and downstream ends of the laboratory channel. The channel had a trapezoidal cross section with 1H:1V side slopes, smooth walls, and a constant bottom width of 1.829 m. The channel dimensions and hydraulic properties are presented in Table 2.
Water 2020, 12, x FOR PEER REVIEW 12 of 76 boundary of the physical model ended up being significantly different than the simulated velocity distribution across this same location. Therefore, the simulated velocity distributions did not begin to match measured velocities until the downstream end of the first channel curve. Simulated and measured velocity distributions matched reasonably well through the second curve. This example illustrates the need for the upstream boundary of a physical model to be well upstream of the region where measurements are of great interest. The details of the 2D and 3D numerical model simulations and the comparison of simulated and measured velocity distributions are presented in Appendix B.

Model Shear Stress Validation
According to verifications by the model developer (Lai 2003a and2003b) and this research, the two-and three-dimensional numerical models (SHR-2D and U2RANS) were considered to be suitable for predicting secondary circulation and the distribution of stream velocity through channel bends. Shear stress computed by the U 2 RANS model was validated for this research by simulating flow through physical models at the University of Iowa [58] and Massachusetts Institute of Technology [59]. The 3D numerical model simulations of shear stress compared well with physical model measurements.

University of Iowa Physical Model
Lai et al. [48] previously verified the U 2 RANS model by comparing simulated water depth and velocity profiles with measurements from a physical model of open channel flow at the University of Iowa [58]. These numerical model simulations were repeated for this research to also verify the model's ability to simulate boundary shear stress. The channel of the physical model was 35.348 m long and included two identical 90° curves of opposite direction (right curve then left curve) connected by a short straight channel segment ( Figure 8). Short, straight channel segments were also included at the upstream and downstream ends of the laboratory channel. The channel had a trapezoidal cross section with 1H:1V side slopes, smooth walls, and a constant bottom width of 1.829 m. The channel dimensions and hydraulic properties are presented in Table 2.   Lai [48] applied the U 2 RANS model to simulate the depth and velocity distributions using the following three different mesh sizes: coarse (45 × 24 × 12), medium (90 × 48 × 24), and fine (180 × 96 × 48). The longitudinal slope of the physical model was mild enough that a flat rigid lid was assumed for the water surface. The flow depths and water surface elevations throughout the model domain were computed from the pressure head simulated by the model at the rigid lid. Super elevations of water depth were measured at cross sections within the physical model. Lai [48] concluded that the numerical simulation results of water depth, using all three mesh sizes, matched the depth measurements very well.
For this study, 49 simulated velocity profiles (using the same coarse, medium, and fine mesh sizes) were compared with measured velocity profiles, at the downstream end of the first channel curve (Figures 9-13). For the simulation using the fine mesh (180 × 96 × 48), the mean dimensionless error between measured and simulated velocities at cross-section S0 was +1.0% (+0.009). The mean error at individual velocity profiles ranged from −10.2% to +6.9% to (−0.091 to +0.071). The overall RMS error between measured and simulated velocities was 6.8% (0.067). The RMS error at individual velocity profiles ranged from 7.3% to 10.6% (0.076 to 0.094). The simulated velocity profiles using all three mesh sizes compared well with the measured profiles. However, the simulated velocities, using the finest model mesh, matched measured profiles best near the channel bottom.        In addition to the simulation of velocity distributions, the U 2 RANS model simulated the boundary shear stress ( Figure 14). Model simulations indicate that boundary shear stress is initially greatest along the right (inside) bank of the first (upstream) channel curve, and then transfers to the left bank of the downstream curve. The maximum shear stress is located along the inside bank near In addition to the simulation of velocity distributions, the U 2 RANS model simulated the boundary shear stress ( Figure 14). Model simulations indicate that boundary shear stress is initially greatest along the right (inside) bank of the first (upstream) channel curve, and then transfers to the left bank of the downstream curve. The maximum shear stress is located along the inside bank near the upstream end of the second channel curve. In addition to the simulation of velocity distributions, the U 2 RANS model simulated the boundary shear stress ( Figure 14). Model simulations indicate that boundary shear stress is initially greatest along the right (inside) bank of the first (upstream) channel curve, and then transfers to the left bank of the downstream curve. The maximum shear stress is located along the inside bank near the upstream end of the second channel curve. Measured and simulated boundary shear stress were compared at four channel cross sections ( Figure 8) through the downstream channel curve (Figures 15-18). For the fine mesh, the overall mean error between measured and simulated shear stress was −2.5% (−0.00008). The mean error at individual cross sections ranged from −1.3% to −4.2% (−0.00005 to −0.00017). For the fine mesh, the overall RMS error between measured and simulated shear stress was 8.0% (0.00032). The RMS error at individual cross sections ranged from 5.2% to 11.2% (0.00021 to 0.00045).
Measured and simulated boundary shear stress were compared at four channel cross sections ( Figure 8) through the downstream channel curve (Figures 15-18). For the fine mesh, the overall mean error between measured and simulated shear stress was −2.5% (−0.00008). The mean error at individual cross sections ranged from −1.3% to −4.2% (−0.00005 to −0.00017). For the fine mesh, the overall RMS error between measured and simulated shear stress was 8.0% (0.00032). The RMS error at individual cross sections ranged from 5.2% to 11.2% (0.00021 to 0.00045). Figure 15. Comparison of simulated and measured shear stress for cross-section CIIO (upstream entrance to second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% percent of maximum velocity).

Figure 15.
Comparison of simulated and measured shear stress for cross-section CIIO (upstream entrance to second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% percent of maximum velocity).
Water 2020, 12, x FOR PEER REVIEW 18 of 76 Figure 16. Comparison of simulated and measured shear stress for cross-section π/16 (11.25° through the second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% of maximum velocity).

Figure 16.
Comparison of simulated and measured shear stress for cross-section π/16 (11.25 • through the second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% of maximum velocity).

Figure 16.
Comparison of simulated and measured shear stress for cross-section π/16 (11.25° through the second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% of maximum velocity).

Figure 17.
Comparison of simulated and measured shear stress for cross-section π/8 (22.5° through the second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% of maximum velocity).

Figure 17.
Comparison of simulated and measured shear stress for cross-section π/8 (22.5 • through the second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% of maximum velocity).
Water 2020, 12, x FOR PEER REVIEW 19 of 76 Figure 18. Comparison of simulated and measured shear stress for cross-section π/4 (45° through the second channel curve). The channel cross section is also plotted along with the location of greatest velocity (within 90% of maximum velocity).
There was no apparent local decrease in the magnitude of measured shear stress at the bottom corners of the trapezoidal cross section. The U 2 RANS model is not able to provide accurate shear stress results at the corners of the trapezoid, therefore, simulated shear stress results were not used from the mesh cells at or next to the bottom corners of the cross section. Otherwise, the simulated shear stress matched measured shear stress quite well. Simulated shear stress was not sensitive to the mesh density along the channel bottom. However, the simulated shear stress did sometimes vary with mesh density along the channel banks. There was no apparent local decrease in the magnitude of measured shear stress at the bottom corners of the trapezoidal cross section. The U 2 RANS model is not able to provide accurate shear stress results at the corners of the trapezoid, therefore, simulated shear stress results were not used from the mesh cells at or next to the bottom corners of the cross section. Otherwise, the simulated shear stress matched measured shear stress quite well. Simulated shear stress was not sensitive to the mesh density along the channel bottom. However, the simulated shear stress did sometimes vary with mesh density along the channel banks.

Massachusetts Institute of Technology Physical Model
An open channel, physical model was developed at the Massachusetts Institute of Technology that had a trapezoidal cross section (2H:1V side slopes) and one 60-degree curve [59]. The width/depth ratio of the laboratory channel was 12 and the Froude number was 0.53. Measured boundary shear stresses from this physical model were reported in the form of a contour plot ( Figure 19). Ippen and Drinker [59] did not provide the individual shear stress measurements that were used to develop the contour plot. The U 2 RANS model was used to simulate the velocity distribution and shear stress of this physical model for a qualitative comparison with measurements ( Figure 19). Again, the longitudinal slope of this physical model slope was mild enough that the slope was neglected, and a flat rigid lid was assumed for the water surface. The patterns of simulated shear stress contours matched rather well with the patterns of the measured contours.

Numerical Model Mesh Development
The necessary length of the straight channel segments, upstream from the meandering reach, was tested by numerical model simulation of an example straight reach from Simulation Set 6 with a length of 80 channel widths ( Figure 20). The simulated bottom shear stress was highest at the upstream boundary (because of the assumed uniform velocity distribution). Then, the shear stress decreased to a minimum after a downstream distance of about two channel widths, and then approached a near constant value after a distance of about 20 channel widths. However, simulation of a longer and separate straight reach was sometimes necessary to determine the reference bottom

Numerical Model Mesh Development
The necessary length of the straight channel segments, upstream from the meandering reach, was tested by numerical model simulation of an example straight reach from Simulation Set 6 with a length of 80 channel widths ( Figure 20). The simulated bottom shear stress was highest at the upstream boundary (because of the assumed uniform velocity distribution). Then, the shear stress decreased to a minimum after a downstream distance of about two channel widths, and then approached a near constant value after a distance of about 20 channel widths. However, simulation of a longer and separate straight reach was sometimes necessary to determine the reference bottom shear stress. The hydraulics through five consecutive meander bends were simulated to ensure repeatability of results. The hydraulic results from the third bend were used to characterize the channel and compared with results from the fourth bend to verify repeatability. The hydraulic conditions through the upstream and downstream most meander bends are influenced to some degree by the straight entrance and exit reaches. Hydraulic results from the middle three of five bends were generally repeatable and not influenced by the upstream or downstream boundary conditions or the straight channel reaches (Figure 21). The experience gained in applying the U 2 RANS model to the CSU physical model was used to guide the 3D model mesh development for the twelve sets of synthetic meandering channels. Four different 3D mesh sizes were used to numerically simulate the physical model. Simulated velocity profiles using the finest mesh (492 × 42 × 32) were nearly the same as results using the next smaller size mesh (327 × 42 × 32) (Appendix B). Therefore, the physical model dimensions (Table 3) were compared with the mesh size from the finest mesh (Table 4) to develop dimensionless relationships.
The hydraulics through five consecutive meander bends were simulated to ensure repeatability of results. The hydraulic results from the third bend were used to characterize the channel and compared with results from the fourth bend to verify repeatability. The hydraulic conditions through the upstream and downstream most meander bends are influenced to some degree by the straight entrance and exit reaches. Hydraulic results from the middle three of five bends were generally repeatable and not influenced by the upstream or downstream boundary conditions or the straight channel reaches ( Figure 21). Water 2020, 12, x FOR PEER REVIEW 22 of 76 The experience gained in applying the U 2 RANS model to the CSU physical model was used to guide the 3D model mesh development for the twelve sets of synthetic meandering channels. Four different 3D mesh sizes were used to numerically simulate the physical model. Simulated velocity profiles using the finest mesh (492 × 42 × 32) were nearly the same as results using the next smaller size mesh (327 × 42 × 32) (Appendix B). Therefore, the physical model dimensions (Table 3) were compared with the mesh size from the finest mesh (Table 4) to develop dimensionless relationships. As previously stated, the U 2 RANS model could not reliably compute shear stress at and near the corner mesh cells of the trapezoidal channels and results from those cells were not used. The   Table 4. Finest numerical model mesh size used to simulate the flow through the CSU physical model (i, j, and k were the model mesh sizes in the longitudinal, transverse, and vertical dimensions, respectively; L, T, and D were the channel dimensions of length, top width, and depth, respectively). As previously stated, the U 2 RANS model could not reliably compute shear stress at and near the corner mesh cells of the trapezoidal channels and results from those cells were not used. The maximum boundary shear stress was found to occur on the channel bottom near the bank toe and, occasionally, on the channel bank above the toe ( Figure 22). Water 2020, 12, x FOR PEER REVIEW 23 of 76 maximum boundary shear stress was found to occur on the channel bottom near the bank toe and, occasionally, on the channel bank above the toe ( Figure 22). A theoretical analysis of the boundary shear stress distribution in straight trapezoidal channels was performed by Olsen and Florey [60] who utilized a membrane analogy. These results indicated that the shear stress would be zero at the bottom corners of a trapezoidal channel. However, data presented by the ASCE Task Committee on Hydraulics, Bank Mechanics, and Modeling of River Width Adjustment [61] indicate that the shear stress is somewhat higher near the bottom corners of a trapezoidal channel and not zero. Ippen and Drinker [59] measured the boundary shear stress in curved trapezoidal channels and found that the maximum shear stress was near the bottom corners of the trapezoidal channel, either along the inside curved bank or along the outside bank at the end of the curve (Figure 19).

Channel Segment
Using a fine mesh size in the model simulations provided more reliable shear stress values closer to the cross-section corners than simulations with a coarse mesh. Therefore, a fine mesh with about one million mesh cells was used to simulate the synthetic meandering channels. The following equations were used to determine the number of mesh cells for the synthetic channels where i, j, and k were the model mesh sizes in the longitudinal, transverse, and vertical dimensions: These equations were based on the ratios of the straight entrance segment (Table 4) with a 12.5% increase to produce about 1 million mesh cells.

Dimensionless Hydraulic and Channel Parameters
Dimensionless parameters were used in the empirical analysis of the 3D model simulation results. Dimensionless model shear stress (τ R ) was computed at each boundary mesh cell as the ratio of the local shear stress at the curved boundary (τ W ) to the equilibrium shear stress (τ o ) along the centerline of a long straight channel with the same cross-sectional shape. This equilibrium shear stress is unique for each simulation as: Dimensionless velocity (V R ) was computed at each cell of model mesh flow field as the ratio of local velocity at a cell to the average cross-sectional channel velocity.
Dimensionless ratios of channel width (W/D) and curvature (W/R c ) were also used in the empirical analysis where W is the wetted channel width, D is the flow depth, and R c is the radius of curvature at the channel centerline. The minimum radius of curvature was used in the empirical correlations. The channel length (l) is made dimensionless by the channel length through one-half of a meander wavelength (Pλ/2), where P is the sinuosity, and λ is the meander wavelength.

Pre-and Post-Processing Programs
Although the 2D and 3D hydraulic models were previously developed, new pre-and post-processing programs were developed for this research. The pre-processing program used Equations (9) and (10) to compute four meandering channel alignments for each simulation (along the top left and top right channel banks and along the left and right channel bank toes). Model meshes were automatically developed from these channel alignments.
The post-processing program computed the following dimensionless hydraulic parameters from the 3D model simulations: In addition, the post-processing program computed the lateral distributions of dimensionless shear stress and identified zones of highest velocity at selected cross sections. Vertical, streamwise velocity profiles were also computed at these selected cross sections at the left bank toe, channel centerline, rank bank toe, and at locations half-way between the channel centerline and bank toes. Although the models can provide results with nine digits of precision, only three digits of precision were used for empirical analysis, consistent with the precision of model input data.

Methodology Limitations
The research methodology was carefully designed to evaluate the increased shear stress through meandering channels. The use of synthetic channels and multidimensional numerical models allowed the variables to be highly controlled. Even though the 2D and 3D models have been verified and validated with physical model measurements, they are not perfect and simulated results do include some error. However, the biggest limitation in the research methodology is the assumption of trapezoidal cross sections with flat bottoms (no slope in the transverse direction). A variable cross-bottom slope could be used to represent the bar-pool channel pattern that creates an additional component of secondary flow [16]. Although the flat-bottom assumption was considered to be a necessary first step in this research, the effects of a variable cross-bottom slope will eventually need to be investigated. Results from this research can help inform hypotheses that can be tested with data from physical model experiments and natural channels, a step that is beyond the scope of this research.

Numerical Model Simulation Results for Meandering Channels
The 2D and 3D hydraulic models were first applied to Simulation Sets 4 through 9 (54 separate model simulations). Model results from these six simulation sets (nine sinuosity in each set) were used to develop empirically based relationships of dimensionless shear stress as a function of the dimensionless channel parameters. Model results from simulation Then, Simulation Sets 1, 2, and 3 and Simulation Sets 10, 11, and 12 were used to validate the empirical equations. Graphs of model results from Simulation Set 5, with a channel sinuosity of 2.75, are presented as a representative example of results produced for all 54 simulations.

Example 2D Model Results
A plan view example of the 2D model mesh is presented in Figure 23 where flow enters at the southwest end. A close-up view of this mesh is presented in Figure 24. Example water-surface elevation contours from the 2D model results are presented in Figure 25.

Example 3D Model Results
The structured 3D model mesh is similar to the 2D mesh, but the vertical mesh lines radiate from along the channel bottom and tend to follow the angle of channel side slopes away from the channel center ( Figure 26). For the example simulation, the structured 3D model mesh has 1587 cells in the

Example 3D Model Results
The structured 3D model mesh is similar to the 2D mesh, but the vertical mesh lines radiate from along the channel bottom and tend to follow the angle of channel side slopes away from the channel center ( Figure 26). For the example simulation, the structured 3D model mesh has 1587 cells in the longitudinal streamwise direction (i-dimension), 31 cells in the cross-stream direction (j-dimension), and 20 cells in the vertical direction (k-dimension) for a total of 983,940 mesh cells.

Example 3D Model Results
The structured 3D model mesh is similar to the 2D mesh, but the vertical mesh lines radiate from along the channel bottom and tend to follow the angle of channel side slopes away from the channel center ( Figure 26). For the example simulation, the structured 3D model mesh has 1587 cells in the longitudinal streamwise direction (i-dimension), 31 cells in the cross-stream direction (j-dimension), and 20 cells in the vertical direction (k-dimension) for a total of 983,940 mesh cells. A color contour shading of surface velocities from the 3D model results is presented in Figure  27. A close-up view of cross-sectional slices of velocity magnitude through a meander bend is presented in Figure 28. A color contour shading of surface velocities from the 3D model results is presented in Figure 27. A close-up view of cross-sectional slices of velocity magnitude through a meander bend is presented in Figure 28.     The thread of maximum velocity alignment is presented in Figure 29. A uniform velocity distribution is specified across the upstream model boundary, and therefore identifying the thread of maximum velocity is difficult through the upstream straight reach. At the beginning of each meander bend, the thread of maximum velocity is near the inside curve and gradually migrates to the channel centerline near the apex (minimum radius of curvature) of the meander bend. Farther downstream the thread of maximum velocity crosses over to the outside of the curve through the downstream portion of the meander bend and remains along that bank until the curvature is reversed and the radius becomes closer to the minimum.
The patten of boundary shear stress contours ( Figure 30) is similar to the pattern of surface velocity contours through the meander bends ( Figure 27).
Dimensionless shear stress, as a function of dimensionless channel length, is presented in Figure 31. The lateral distribution of simulated shear stress and maximum velocity were compiled for four selected cross sections and their locations are presented in Figure 32. Cross-section plots of shear stress and zones of maximum velocity for cross-sections 112, 702, 810, and 950 are presented in Figures 33-36. Vertical velocity profiles plots are presented for these four cross sections are presented in  downstream the thread of maximum velocity crosses over to the outside of the curve through the downstream portion of the meander bend and remains along that bank until the curvature is reversed and the radius becomes closer to the minimum. The patten of boundary shear stress contours ( Figure 30) is similar to the pattern of surface velocity contours through the meander bends ( Figure 27).
Dimensionless shear stress, as a function of dimensionless channel length, is presented in Figure  31. The lateral distribution of simulated shear stress and maximum velocity were compiled for four selected cross sections and their locations are presented in Figure 32 Model results of velocity and shear stress were different through the first meander bend, but results were repeatable through subsequent meander bends. Therefore, model results were not used from the first meander bend, but from the second-to-last meander bend in the series. The lateral distributions of velocity and boundary shear stress were symmetrical about the channel centerline through the upstream straight reach. These symmetrical distributions caused the maximum nearbank shear stress through the first meander bend to be lower than through subsequent meander bends.
By the downstream end of the first meander bend, velocity and shear stress increased along the left (outside) bank of the first bend, which is the inside bank at the beginning of the second meander bend. The zone of high velocity and shear stress continued along the inside bank of the next meander bend and reached a maximum before the bend apex where the radius of curvature was at a minimum.                     The longitudinal distance required for the thread of high velocity and shear stress to be transferred from the inside bank to the outside bank is known as the planform phase lag, which is found in natural channels. The simulated phase lags were long enough that the zone of highest velocity and shear stress was along the outside bank near the downstream end of the meander bend, which is also along the inside bank at the beginning of the next (downstream) meander bend. The laboratory experiments by Yen [58] and Ippen and Drinker [59] also demonstrated that the zone of highest shear stress was along the inside channel bank through the curve, and then transferred to the outside bank near the downstream end of the curve.
The plan view and cross-section plots (Figures 32-36) show that boundary shear stress near the Model results of velocity and shear stress were different through the first meander bend, but results were repeatable through subsequent meander bends. Therefore, model results were not used from the first meander bend, but from the second-to-last meander bend in the series. The lateral distributions of velocity and boundary shear stress were symmetrical about the channel centerline through the upstream straight reach. These symmetrical distributions caused the maximum near-bank shear stress through the first meander bend to be lower than through subsequent meander bends.
By the downstream end of the first meander bend, velocity and shear stress increased along the left (outside) bank of the first bend, which is the inside bank at the beginning of the second meander bend. The zone of high velocity and shear stress continued along the inside bank of the next meander bend and reached a maximum before the bend apex where the radius of curvature was at a minimum. By the downstream end of the second meander bend (curving to the left), secondary currents began to increase velocity and shear stress along the opposite, right bank. The transfer of high velocity and shear stress from the inside bank to the outside bank was repeated through subsequent meander bends (Figures 27, 29 and 30).
The longitudinal distance required for the thread of high velocity and shear stress to be transferred from the inside bank to the outside bank is known as the planform phase lag, which is found in natural channels. The simulated phase lags were long enough that the zone of highest velocity and shear stress was along the outside bank near the downstream end of the meander bend, which is also along the inside bank at the beginning of the next (downstream) meander bend. The laboratory experiments by Yen [58] and Ippen and Drinker [59] also demonstrated that the zone of highest shear stress was along the inside channel bank through the curve, and then transferred to the outside bank near the downstream end of the curve.
The plan view and cross-section plots (Figures 32-36) show that boundary shear stress near the banks increases due to channel curvature. At the channel crossing location, where the zone of high velocity and shear stress is transferring from the inside bank to the outside bank, the shear stress is somewhat uniform across the channel bottom, but the velocity distribution is not uniform because of the secondary currents. The vertical velocity profiles show the model is capable of computing the maximum velocity, at a given profile, that is often at the water surface, but sometimes below the water surface (Figures 37-40).
The zones of high velocity and shear stress, alternating along the right and left channel banks, would tend to erode the channel along those locations and from pools with point bars depositing along the opposite, inside bank. This would result in a transverse channel bed slope. However, the model simulations have fixed bed and bank boundaries and erosion is not simulated. Seminara [16] pointed out that for channels with non-erodible banks, flow velocity tended to be higher along the inner bank. Farther downstream, the secondary currents transfer momentum toward the outer bend, moving the thread of high velocity from the inner bank to the outer bank.
The model results from all simulation sets indicate a similar longitudinal distribution of shear stress. The longitudinal increase in shear stress along a channel bank is largely parabolic. The longitudinal locations where the shear stress along both banks is the same corresponds closely with the beginning and ending of each meander curve. This is illustrated by the examples from Simulation Sets 7, 8, and 9 for a sinuosity of 2.00 (Figures 41-43). Model results indicate that the location of maximum shear stress is one-quarter to one-half of the way through the next meander bend, but upstream of the bend apex where the radius of curvature reaches a minimum. The maximum shear stress tends to reduce with an increase in the width/depth ratio. maximum shear stress is one-quarter to one-half of the way through the next meander bend, but upstream of the bend apex where the radius of curvature reaches a minimum. The maximum shear stress tends to reduce with an increase in the width/depth ratio.   maximum shear stress is one-quarter to one-half of the way through the next meander bend, but upstream of the bend apex where the radius of curvature reaches a minimum. The maximum shear stress tends to reduce with an increase in the width/depth ratio.

Shear Stress Relationships
Boundary shear stress results from 3D model simulations of Simulation Sets 4 through 9 were compiled and evaluated to determine empirical relationships with channel curvature, width, and depth. The phase lag location and magnitude of maximum and minimum shear stress due to channel curvature was determined for each simulation (Figures 41-43). The shear stress results from the curvature of the third consecutive meander bend were used to avoid any influence from the upstream and downstream model boundaries.

Shear Stress Magnitude
Maximum dimensionless shear stresses, along the left and right channel banks for Simulation Sets 4 through 9, varied between 1.21 and 1.54, while minimum shear stress varied between 0.46 and 0.86. Maximum shear stresses were highest, and minimum sinuosity were lowest, at sinuosity between 1.50 and 2.25 ( Figure 44). However, shear stresses did not correlate well with sinuosity.
Correlations were found between the maximum dimensionless near-bank shear stress and the ratio of channel top to minimum radius of curvature ( Figure 45). The correlations were found by segregating the data by the channel width/depth ratio (W/D) as follows: All three linear regression lines for these width/depth ratio sets had similar positive slopes and correlation coefficients (R 2 ) of 0.94, 0.60, and 0.68. Maximum shear stress tended to decrease with increases in the width/depth ratio.

Shear Stress Relationships
Boundary shear stress results from 3D model simulations of Simulation Sets 4 through 9 were compiled and evaluated to determine empirical relationships with channel curvature, width, and depth. The phase lag location and magnitude of maximum and minimum shear stress due to channel curvature was determined for each simulation (Figures 41-43). The shear stress results from the curvature of the third consecutive meander bend were used to avoid any influence from the upstream and downstream model boundaries.

Shear Stress Magnitude
Maximum dimensionless shear stresses, along the left and right channel banks for Simulation Sets 4 through 9, varied between 1.21 and 1.54, while minimum shear stress varied between 0.46 and 0.86. Maximum shear stresses were highest, and minimum sinuosity were lowest, at sinuosity between 1.50 and 2.25 ( Figure 44). However, shear stresses did not correlate well with sinuosity.
Correlations were found between the maximum dimensionless near-bank shear stress and the ratio of channel top to minimum radius of curvature ( Figure 45). The correlations were found by segregating the data by the channel width/depth ratio (W/D) as follows:       All three linear regression lines for these width/depth ratio sets had similar positive slopes and correlation coefficients (R 2 ) of 0.94, 0.60, and 0.68. Maximum shear stress tended to decrease with increases in the width/depth ratio.
The best correlation (r 2 = 0.94) was found for narrow channels with width/depth ratios between 4 and 10: Therefore, this initial regression equation was used to compute the maximum shear stress for all 54 channels. The dimensionless residual between the maximum shear stress, computed by Equation (17), and the maximum shear stress simulated by the 3D model were plotted against the channel width/depth ratio ( Figure 46). For the narrow channels (width/depth ratios between 4 and 10), the dimensionless residuals were low (−0.015 to +0.030). For the wider channels, the data were segregated into three ranges of width/depth ratio and a series of linear relationships were found between the residuals and the width/depth ratios. For width/depth ratios between 10 and 40, the residuals decreased (in the negative direction) with increases in width/depth ratio. For width/depth ratios greater than 40, residuals did not corelate with the width/depth ratio. The best correlation (r 2 = 0.94) was found for narrow channels with width/depth ratios between 4 and 10: Therefore, this initial regression equation was used to compute the maximum shear stress for all 54 channels. The dimensionless residual between the maximum shear stress, computed by Equation (17), and the maximum shear stress simulated by the 3D model were plotted against the channel width/depth ratio ( Figure 46). For the narrow channels (width/depth ratios between 4 and 10), the dimensionless residuals were low (−0.015 to +0.030). For the wider channels, the data were segregated into three ranges of width/depth ratio and a series of linear relationships were found between the residuals and the width/depth ratios. For width/depth ratios between 10 and 40, the residuals decreased (in the negative direction) with increases in width/depth ratio. For width/depth ratios greater than 40, residuals did not corelate with the width/depth ratio.  (17), and the maximum shear stress simulated by the 3D model versus the channel width/depth ratio.
The regression equations of residuals were combined with Equation (17) to produce the following equation to estimate the maximum shear stress through a meandering trapezoidal channel:

Shear Stress Phase Lag
In all 3D model simulations, the phase lag location of maximum shear stress occurred downstream from the end of the meander bend and through the beginning of the next meander bend where the channel curvature reverses direction. The shear stress phase lag is defined as the channel length between the point of maximum curvature (minimum radius at the bend apex) and the downstream location of maximum near-bank shear stress. The simulated dimensionless phase lags from the 3D ranged from 0.845 to 1.075. The dimensionless shear phase lag was found to correlate with the channel sinuosity ( Figure 47). The phase lag data were also segregated by the channel width/depth ratio (W/D) into the following three categories:

Shear Stress Phase Lag
In all 3D model simulations, the phase lag location of maximum shear stress occurred downstream from the end of the meander bend and through the beginning of the next meander bend where the channel curvature reverses direction. The shear stress phase lag is defined as the channel length between the point of maximum curvature (minimum radius at the bend apex) and the downstream location of maximum near-bank shear stress. The simulated dimensionless phase lags from the 3D ranged from 0.845 to 1.075. The dimensionless shear phase lag was found to correlate with the channel sinuosity ( Figure 47). The phase lag data were also segregated by the channel width/depth ratio (W/D) into the following three categories:  Most 3D model simulation results show a trend of decreasing phase lag with increases in channel sinuosity. However, phase lags tended to increase for sinuosity values greater than 2.25 for narrow channels with a width/depth ratio less than 10. For a given sinuosity, the phase lag was longest for channels with width/depth ratios less than 10. For width/depth ratios between 4 and 10, a secondorder polynomial regression equation was found to have a correlation coefficient (R 2 ) of 0.87. This equation predicts the phase lag is longest for a sinuosity of 1.1 and at a minimum for a sinuosity of 2.25.  Most 3D model simulation results show a trend of decreasing phase lag with increases in channel sinuosity. However, phase lags tended to increase for sinuosity values greater than 2.25 for narrow channels with a width/depth ratio less than 10. For a given sinuosity, the phase lag was longest for channels with width/depth ratios less than 10. For width/depth ratios between 4 and 10, a second-order polynomial regression equation was found to have a correlation coefficient (R 2 ) of 0.87. This equation predicts the phase lag is longest for a sinuosity of 1.1 and at a minimum for a sinuosity of 2.25.
For width depth ratios between 10 and 20, a linear regression equation was found to have a correlation coefficient (R 2 ) of 0.89. This equation predicts that the phase lag is longest for a sinuosity of 1.10 and decreases with increases in sinuosity.
Phase Lag = −0.0763P + 1.070, 10 < W /D ≤ 20 (20) For width/depth ratios between 20 and 84, the phase lag is nearly constant with sinuosity and the average phase lag is 0.908. Because the phase lag only varies between 0.878 and 0.942 (between −3.3% and +3.7% of the mean), the slope of the linear regression equation is relatively flat, so the correlation coefficient (R 2 ) is only 0.32: The maximum shear stress and phase lag results from the 3D numerical modeling and predicted values from the empirical equations are presented in Appendix C.

Validation and Uncertainty
The empirical regressions equations (Equations (18)-(21)) were validated by application to the six simulation sets with the lowest and greatest discharge (Simulation Sets 1, 2, 3, and 10, 11,12) and meandering channels with the steepest and most mild slopes (Simulation Sets 3 and 10). Then, the predictions from the empirical equations were compared with the 3D model results (Tables 5 and 6) from the six simulation sets not used in the regression analyses. Each of these six simulation sets included three river channels with sinuosity values of 1.75, 2.25, and 2.75 for a total of 18 synthetic meandering river channels. Table 5. Validation of empirical equation to predict maximum shear stress for Simulation Sets 1-3 and 10-12. The differences in shear stress (simulated by the 3D model and computed from Equation (18)) were computed for each channel and the minimum and maximum differences are reported here. Empirical predictions of maximum dimensionless shear stress and phase lag agreed well with the 3D model results not used in the regression analyses. Empirical predictions were within 1% of 3D model results for shear stress and within 3% for phase lag. RMS errors were 0.021 (1.4% of average) and 0.027 (2.9% of average).

Maximum Dimensionless Shear
Given that the empirical equations agreed so well with 3D model results not used in the regression, new regression equations were developed using simulations from all twelve sets (a total of 72 synthetic meandering river channels). The same regression approach was used as before and coefficients for the new equations were close in magnitude to the coefficients of the previous equations.

Shear Stress Magnitude
Using narrow river channels with width/depth ratios between 4 and 10 (from Simulation sets 1, 4, 7, 9, and 12), linear regression produced the following equation to predict the maximum shear stress ( Figure 48): Water 2020, 12, x FOR PEER REVIEW 42 of 76 Given that the empirical equations agreed so well with 3D model results not used in the regression, new regression equations were developed using simulations from all twelve sets (a total of 72 synthetic meandering river channels). The same regression approach was used as before and coefficients for the new equations were close in magnitude to the coefficients of the previous equations.

Shear Stress Magnitude
Using narrow river channels with width/depth ratios between 4 and 10 (from Simulation sets 1, 4, 7, 9, and 12), linear regression produced the following equation to predict the maximum shear stress ( Figure 48): This equation was used to compute the maximum shear stress for all 72 simulations. The residuals between the maximum shear stress, computed by Equation (22), and the maximum shear stress simulated by the 3D model were correlated with the width/depth ratio ( Figure 49).  This equation was used to compute the maximum shear stress for all 72 simulations. The residuals between the maximum shear stress, computed by Equation (22), and the maximum shear stress simulated by the 3D model were correlated with the width/depth ratio ( Figure 49).  Then, the regression equations of residuals were combined with Equation (22), to produce the following new equation to estimate the maximum shear stress through a meandering trapezoidal channel: Equation (23) was applied to the data from all twelve simulation sets to estimate the maximum dimensionless shear stress. The empirical estimates of maximum shear stress were correlated with results from the 3D model simulations to determine the 5 and 95% confidence limits about the estimated values (±0.031) ( Figure 50).  (23)) and simulated by the 3D model.

Validation of Shear Stress Phase Lag
The phase lag of the modeled shear stress, from all twelve simulation sets, was correlated with channel sinuosity (Figure 51). The new linear regression lines (solid lines) match very closely with the initial regression lines (dashed lines) that were developed using data from half the simulation sets.  (23)) and simulated by the 3D model.

Validation of Shear Stress Phase Lag
The phase lag of the modeled shear stress, from all twelve simulation sets, was correlated with channel sinuosity (Figure 51). The new linear regression lines (solid lines) match very closely with the initial regression lines (dashed lines) that were developed using data from half the simulation sets. Figure 50. Comparison of maximum dimensionless near-bank shear stresses (computed by Equation (23)) and simulated by the 3D model.

Validation of Shear Stress Phase Lag
The phase lag of the modeled shear stress, from all twelve simulation sets, was correlated with channel sinuosity (Figure 51). The new linear regression lines (solid lines) match very closely with the initial regression lines (dashed lines) that were developed using data from half the simulation sets.   Simulated phase lags ranged from 0.841 to 1.075. These data were again segregated by the channel width/depth ratio (W/D) into the following three categories: New regression equations were developed using 3D model results from all simulation sets. For width/depth ratios between 5 and 10, a second-order polynomial regression equation was found to have a correlation coefficient (R2) of 0.68. This equation predicts that the phase lag is longest for a sinuosity of 1.10 and decreases to a minimum for a sinuosity of 2.25.
For width depth ratios between 10 and 20, a linear regression equation was found to have a correlation coefficient (R 2 ) of 0.88.
Phase Lag = −0.0765P + 1.071, 10 < W /D ≤ 20 (25) For width depth ratios between 20 and 84, the phase lag is nearly constant with sinuosity and the average phase lag is 0.907. Because the phase lag only varied between 0.878 and 0.942 (between −3.2% and +3.9% of the mean), the slope of the linear regression equation was relatively flat, so the correlation coefficient (R 2 ) was only 0.26: A combination of Equations (24)-(26) were used to estimate the shear stress phase lag for each meandering channel of the twelve simulation sets and the results were compared with the phase lags simulated by the 3D model ( Figure 52). The regression slope of the comparison (0.820) was somewhat lower than the line of perfect agreement but had a reasonable correlation coefficient (R 2 = 0.822). The confidence interval about the estimated values is ±0.037, which represents the 5 and 95 percent confidence limits. The maximum shear stress and phase lag results from the 3D numerical modeling and predicted values from the empirical equations are presented in Appendix C. ℎ = −0.0765 + 1.071, 10 < ≤ 20 ( For width depth ratios between 20 and 84, the phase lag is nearly constant with sinuosity the average phase lag is 0.907. Because the phase lag only varied between 0.878 and 0.942 (betw −3.2% and +3.9% of the mean), the slope of the linear regression equation was relatively flat, so correlation coefficient (R 2 ) was only 0.26: ℎ = −0.0140 + 0.936, 20 < ≤ 84 ( A combination of Equations (24)- (26) were used to estimate the shear stress phase lag for e meandering channel of the twelve simulation sets and the results were compared with the phase simulated by the 3D model ( Figure 52). The regression slope of the comparison (0.820) was somew lower than the line of perfect agreement but had a reasonable correlation coefficient (R 2 = 0.822). confidence interval about the estimated values is ±0.037, which represents the 5 and 95 perc confidence limits. The maximum shear stress and phase lag results from the 3D numerical mode and predicted values from the empirical equations are presented in Appendix C.

Discussion
The simulation of a wide range of meandering channels using a 3D numerical model provided consistent datasets for testing the hypothesis that the maximum shear stress could be related to channel curvature and width. The simulation of long straight entrance and exit channels and consecutive meander bends in between was important to avoid the influence of model boundary conditions. Multiple consecutive meander bends demonstrated that simulated results were repeatable. The meandering channel alignments, defined by sine-generated curves, had smooth variations in the radius of curvature that ranged from infinite at the straight channel segments and crossings to minimums at the meander bend apexes. The smooth variations in radius of curvature were more like natural channels than curves of a constant radius [53].
The numerical simulation of river hydraulics produced datasets that represented a wide range of meandering river channels of the type (width and depth) found in nature. The model simulation sets spanned three orders of magnitude in bankfull discharge (2.82, 28.3, 283, and 2830 m 3 /s) and spanned four orders of magnitude in longitudinal slope (2.59 × 10 −6 to 1.54 × 10 −2 ). Channel width/depth ratios ranged from 4.7 to 84 and median sediment grain sized ranged from 0.125 to 4 mm. The use of synthetic channels and simulations with multidimensional numerical models allowed the variables to be highly controlled. The evaluation of dimensionless values (e.g., velocities and shear stresses simulated by the 3D model) reduced the effects of uncertainty in model inputs. For example, increases or decreases in channel roughness would change the shear stress magnitudes, but the change in the dimensionless values was much less.
The case can be made that validated multidimensional numerical models are uniquely suited to produce hydraulic datasets over a wide range of channel slopes and discharges. Shear stress measurements from a wide range of natural channels at near bankfull discharge would be logistically difficult to obtain using consistent methods and each meander bend alignment could have unique characteristics. Channel and discharge characteristics can be controlled in physical models, but the construction of many tens of models with consecutive meander bends and experimental measurements could be rather expensive.
Empirical analysis of the 3D model results indicated that the maximum dimensionless near-bank shear stress is linearly correlated with the ratio of channel width to minimum radius of curvature and the channel width/depth ratio. Maximum shear stress was found to linearly increase with increases in the ratio of width to radius of curvature and decrease with increases in the width/depth ratio. These findings are consistent with the literature [16,17,22] In addition, empirical analysis indicated that the maximum shear stress was downstream from the end of the meander bend and within the next meander bend where there was a reversal of channel curvature. This is consistent with findings by Seminara [16] for channels with non-erodible banks. Long phase lags were also reported for physical models with fixed trapezoidal cross sections [58,59]. The flat bottoms of the trapezoidal channels did not have thalweg pools along the outside of meander bends where discharge and shear stress could become more concentrated [16]. The 3D model phase lag between the location of minimum radius of curvature and maximum shear stress was found to correlate with channel sinuosity and the width/depth ratio. For wide channels (W/D > 10), the phase lag linearly decreased with increases in sinuosity. For wider channels (W/D > 20), the phase lag was nearly constant (0.93). For narrow channels (W/D < 10), the phase lag parabolically varied with sinuosity, with minimum phase lag near a sinuosity of 2.25.
The empirical equations to predict dimensionless near-bank shear stress magnitude and phase lag (Equations (24)- (26)) are limited to the idealized geometry of the meandering river channels with trapezoidal cross sections and alignments of a sign-generated curve. The constant width, trapezoidal cross sections all had 2V:1H side slopes with a flat bottom, which was in contrast to natural meandering channels that typically have a thalweg along the outside of the meander bend and some variation in width [14,16,[25][26][27][28].
The simulation of erodible channels or channels with a variable transverse bottom slope was beyond the scope of this research but would be a good topic for additional investigation. A constructed meandering channel, with a constant trapezoidal cross section and erodible boundaries, would be expected to have high zones of shear stress in response to the meander bends. Over time, bed and bank erosion would be expected to create thalweg pools, point bars, and riffles. Depending on sediment supply and streambank erodibility, channel width could vary with longitudinal distance [14,16,[25][26][27][28]. The evolution of channel geometry through time would be expected to affect the distribution and magnitude of velocity and shear stress [30][31][32][33][34][35][36][37][38]. For example, the formation of point bars and pools would be expected to concentrate discharge and increase velocity and shear stress along the outside of meander bends [16]. Therefore, new empirical equations are needed to account for the variable transverse bottom slope of the channel cross section.
A future extension of this research would be to expand a few of the simulation sets reported here (e.g., Simulation sets 7, 8, and 9) to include channels with a variable transverse slope to simulate pools, bars, and crossings associated the meandering channels. The transverse slope of the channel bottom could gradually increase in the downstream direction from zero at the channel crossing to a maximum transverse slope at some location downstream from the bend apex (minimum radius of curvature). From there, the transverse slope could gradually decrease to the next downstream crossing where the transverse slope would reverse direction through the next meander bend. Physical model experiments with erodible channel bottoms would help guide the development of synthetic channels for 3D numerical modeling. A range of maximum transverse slopes could be simulated as a function of sediment grain size. For each synthetic channel, the locations of maximum transverse slope (deepest pool) would have to be selected iteratively so there is a reasonable match with the resulting phase lag simulated by the 3D model. The channel crossing locations could be set at halfway between the deepest pool locations.

Conclusions
The application of a verified multidimensional numerical hydraulics model to develop detailed datasets, representing a wide range of river channels, is a useful method for studying hydraulic relationships. This method produced quantitative relationships between the maximum dimensionless shear stress, and the phase lag location, with the dimensionless channel ratios of curvature, width, and sinuosity. The increase in shear stress due to curvature was found to be independent of discharge, slope, and grain size.
Numerical modeling does not replace the need for detailed hydraulic measurements from physical model experiments or natural channels. Laboratory experiments and field measurements, along with multidimensional numerical modeling, can be quite complementary. The use of multidimensional models can provide greater efficiency than only physical model experiments or analyses of field data measurements.
Funding: This research received no external funding.
Acknowledgments: James C.Y. Guo (University of Colorado at Denver) provided advice and encouragement throughout the research. Yong G. Lai (Bureau of Reclamation) graciously provided training and mentoring on the proper application of his 2D and 3D hydraulic models.

Conflicts of Interest:
The author declares no conflict of interest.

2D
Two-dimensional 3D Three-dimensional SRH-2D Two-dimensional numerical model to simulate depth-averaged hydraulics of open channels U2RANS Three-dimensional numerical model to simulate channel hydraulics

Appendix A. Synthetic Meandering Channels
This appendix presents the dimensions and hydraulic properties of the synthetic meandering channels. Channel dimensions are presented in Table A1 while hydraulic properties are presented in Table A2.

Appendix B. Numerical Model Simulations of Physical Model Experiment
This appendix describes the 2D and 3D numerical model simulations of velocity distributions through a large-scale physical model (laboratory flume) at Colorado State University (CSU) [7]. The physical model channel was 57.5 m long and included two channel curves of different constant radii and width, separated by a short transitional straight reach, 6.355 m long [56,57] (Table A3). The physical model had a trapezoidal cross section with 3H:1V side slopes. The upstream channel curve was 125 degrees to the right while the downstream curve was 73 degrees to the left. The channel width and radius of curvature were constant through each channel curve. The channel width of the downstream curve was about three-fourths the width of the upstream curve. The hydraulic characteristics of the physical model are presented in Table A4.  [56,57]. A total of 126 velocity profiles were measured in the physical model ( Figure A1). Measurement cross sections were numbered from 1, at the upstream end, to 18, at the downstream end. Velocity profiles were measured at the following seven locations along each cross section and those measurement locations are labeled "a" through "g", left to right looking downstream: • Measurement locations "a" and "b" were on the left bank side slope; • Measurement location "c" was at the toe of the left bank; • Measurement location "d" was at the center of the channel; • Measurement location "e" was at the toe of the right bank; • Cross-section measurement locations "f" and "g" were on the right bank side slope.
• Measurement locations "a" and "b" were on the left bank side slope; • Measurement location "c" was at the toe of the left bank; • Measurement location "d" was at the center of the channel; • Measurement location "e" was at the toe of the right bank; • Cross-section measurement locations "f" and "g" were on the right bank side slope.

Appendix B.1. Numerical Model Simulation Approach
A structured mesh was developed for the SRH-2D model that included 326 cells in the longitudinal direction and 28 cells in the cross-stream direction for a total of 9128 mesh cells ( Figure A2). Of the 28 mesh cells in the cross-stream direction, 20 cells represented the channel bottom while four cells represented the left and the right channel side slopes ( Figure A3). The numerical model domain of the channel was extended the equivalent of ten channel widths upstream and downstream (at the same slope) past the boundaries of the physical model to prevent the numerical model boundaries from influencing the solution domain representing the physical model. The upstream boundary for SRH-2D was specified as a constant discharge with the assumption of uniform velocity distribution parallel to the channel. The downstream boundary was specified as a level water surface elevation corresponding to normal depth.
(at the same slope) past the boundaries of the physical model to prevent the numerical model boundaries from influencing the solution domain representing the physical model. The upstream boundary for SRH-2D was specified as a constant discharge with the assumption of uniform velocity distribution parallel to the channel. The downstream boundary was specified as a level water surface elevation corresponding to normal depth.      The SRH-2D numerical model was calibrated by adjusting the water surface elevation at the downstream boundary, and the channel roughness, and therefore that the simulated water surface profile matched the measurements from the physical model ( Figure A4). The calibrated, downstream-boundary, water-surface elevation was 29.868 m, which corresponds to a flow depth of 0.320 m. The Manning's, n, roughness coefficient was calibrated to 0.0186, which is consistent with the smooth concrete surface of the physical model. The SRH-2D numerical model was calibrated by adjusting the water surface elevation at the downstream boundary, and the channel roughness, and therefore that the simulated water surface profile matched the measurements from the physical model ( Figure A4). The calibrated, downstreamboundary, water-surface elevation was 29.868 m, which corresponds to a flow depth of 0.320 m. The Manning's, n, roughness coefficient was calibrated to 0.0186, which is consistent with the smooth concrete surface of the physical model. Similar to the 2D model mesh, the 3D model mesh was extended ten channel widths farther upstream and downstream, and at the same longitudinal slope, past the boundaries of the physical model ( Figure A5). The top surface of the 3D mesh corresponds to the simulated water surface from the 2D model. A uniform velocity distribution (laterally and vertically) was specified at the upstream model boundary. Similar to the 2D model mesh, the 3D model mesh was extended ten channel widths farther upstream and downstream, and at the same longitudinal slope, past the boundaries of the physical model ( Figure A5). The top surface of the 3D mesh corresponds to the simulated water surface from the 2D model. A uniform velocity distribution (laterally and vertically) was specified at the upstream model boundary. The SRH-2D and U 2 RNAS models are both robust and compute reasonable results given correct input data. However, the velocity distributions simulated by the U 2 RANS model can be sensitive to the mesh density. Therefore, the sensitivity of the mesh density was evaluated by simulating velocities using four different structured mesh sizes. At least one mesh dimension was kept the same between two consecutive meshes. The number of mesh points along the downstream (I) direction, cross stream (J) direction, and vertical (K) direction are listed below for each of the following four mesh sizes: A close-up plan view of the 492 × 42 × 32 mesh size is presented in Figure A6 and a cross-section view is presented in Figure A7. The cross-sectional mesh lines in the K-direction are vertical at the channel center, and then gradually transition to match the channel bank slopes on the left and rightsides. The SRH-2D and U 2 RNAS models are both robust and compute reasonable results given correct input data. However, the velocity distributions simulated by the U 2 RANS model can be sensitive to the mesh density. Therefore, the sensitivity of the mesh density was evaluated by simulating velocities using four different structured mesh sizes. At least one mesh dimension was kept the same between two consecutive meshes. The number of mesh points along the downstream (I) direction, cross stream (J) direction, and vertical (K) direction are listed below for each of the following four mesh sizes: A close-up plan view of the 492 × 42 × 32 mesh size is presented in Figure A6 and a cross-section view is presented in Figure A7. The cross-sectional mesh lines in the K-direction are vertical at the channel center, and then gradually transition to match the channel bank slopes on the left and right-sides.

Appendix B.2. Comparison of 3D model Velocities with Measurements
For this laboratory experiment, the flow velocities entering the channel at Cross-section 1 ( Figure  A1) were not uniformly distributed and highest along the inside of the channel curve (right side). About halfway through the channel curve (at Cross-section 4), the measured velocities were more uniform. By the downstream end of channel curve (Cross-section 8), measured velocities were greatest along the outside of the curve (left side). By the downstream end of straight transition reach (Cross-section 10), flow velocities were still greatest along the left side of the channel. By the downstream end of the second channel curve (Cross-section 18), the highest flow velocities had transferred to the outside portion of the curve (right side).
Simulated and measured velocities were compared graphically, both in plan and profile views. Each of the 126 measured vertical velocity profiles included five to nine separate velocity measurements (U, V, W). Simulated velocity vectors were interpolated from the 3D model results at the same locations as the measured velocity vectors. An overview plan-view comparison of simulated and measured surface velocity vectors (using the finest 3D mesh) is presented in Figure A8. Additional close-up plan-views of the velocity vector comparisons are presented in Figures A9-A12.

Appendix B.2. Comparison of 3D model Velocities with Measurements
For this laboratory experiment, the flow velocities entering the channel at Cross-section 1 ( Figure  A1) were not uniformly distributed and highest along the inside of the channel curve (right side). About halfway through the channel curve (at Cross-section 4), the measured velocities were more uniform. By the downstream end of channel curve (Cross-section 8), measured velocities were greatest along the outside of the curve (left side). By the downstream end of straight transition reach (Cross-section 10), flow velocities were still greatest along the left side of the channel. By the downstream end of the second channel curve (Cross-section 18), the highest flow velocities had transferred to the outside portion of the curve (right side).
Simulated and measured velocities were compared graphically, both in plan and profile views. Each of the 126 measured vertical velocity profiles included five to nine separate velocity measurements (U, V, W). Simulated velocity vectors were interpolated from the 3D model results at the same locations as the measured velocity vectors. An overview plan-view comparison of simulated and measured surface velocity vectors (using the finest 3D mesh) is presented in Figure A8. Additional close-up plan-views of the velocity vector comparisons are presented in Figures A9-A12.

Appendix B.2. Comparison of 3D Model Velocities with Measurements
For this laboratory experiment, the flow velocities entering the channel at Cross-section 1 ( Figure A1) were not uniformly distributed and highest along the inside of the channel curve (right side). About halfway through the channel curve (at Cross-section 4), the measured velocities were more uniform. By the downstream end of channel curve (Cross-section 8), measured velocities were greatest along the outside of the curve (left side). By the downstream end of straight transition reach (Cross-section 10), flow velocities were still greatest along the left side of the channel. By the downstream end of the second channel curve (Cross-section 18), the highest flow velocities had transferred to the outside portion of the curve (right side).
Simulated and measured velocities were compared graphically, both in plan and profile views. Each of the 126 measured vertical velocity profiles included five to nine separate velocity measurements (U, V, W). Simulated velocity vectors were interpolated from the 3D model results at the same locations as the measured velocity vectors. An overview plan-view comparison of simulated and measured surface velocity vectors (using the finest 3D mesh) is presented in Figure A8         Another way to compare simulated and measured velocities is by plotting the relative depthaveraged, stream-wise velocities versus the measurement cross sections. Velocities at three locations were plotted for each cross section representing the channel center (location d) and the toe of each channel bank (locations c and e) ( Figure A13).
There were both differences and similarities between simulated and measured depth-averaged Another way to compare simulated and measured velocities is by plotting the relative depth-averaged, stream-wise velocities versus the measurement cross sections. Velocities at three locations were plotted for each cross section representing the channel center (location d) and the toe of each channel bank (locations c and e) ( Figure A13).
There were both differences and similarities between simulated and measured depth-averaged velocities due to different velocity distributions at the upstream boundaries of the numerical and physical models. The measured stream-wise velocities at Cross-section 1 were greater along both channel banks than the middle of the channel ( Figure A9). However, the straight channel reach of the numerical model resulted in the simulated velocities being greatest along the channel center at the beginning of the first channel curve.
For the first four cross sections, measured velocities were greater along the right side (inside of curve, measurement location e) than along the left side (outside of curve, measurement location c). These measured velocities along left and right sides of the channel were greater than the measured velocities along the channel centerline (measurement location d). However, by Cross-section 5, the measured velocities were of nearly the same magnitude along the left side (outside of curve), right side (inside of curve), and channel centerline. Through the next two cross sections of the upstream curve (Cross-sections 6 and 7), the highest velocities were along the channel centerline and lowest velocities were along the right side of the channel curve (inside of curve). By the downstream end of this first curve (Cross-section 8), the highest velocity was along the left side of the curve (outside of curve, measurement location c) and the lowest along the right side of the channel curve (inside of curve, measurement location e). This trend of measured velocities continued through the straight channel reach and through the first four cross sections of the next channel curve (Cross-sections 10-13). Measured velocities were greatest along the left side (inside of curve) through Cross-section 12. Left-side velocities began a decreasing trend downstream to Cross-section 17. Measured velocities along the right side (outside of curve) generally increased throughout the length of downstream curve and were greater than velocities along the inside of the curve by Cross-section 16.  In contrast, simulated velocities were greatest along the channel centerline through both channel curves (location d). At the upstream end of the first channel curve (Cross-section 1), simulated velocities were nearly the same along both the left and right sides of the channel curve. As flow continued through the channel curve, simulated velocities were greater along the left side (outside of curve, location c) than along the right side (inside of curve, location e). This trend in simulated velocities continued through the straight reach and the beginning of the downstream curve. However, the reversal in channel curvature of the downstream curve caused the simulated velocities along the right bank (outside of curve) to increase through the entire curve, while the simulated velocities along the left side (inside of curve) decreased.
Because of different velocity distributions at the upstream boundaries of the numerical model and physical model, simulated velocities along the left bank (outside of curve, location c) did not match measurements until Cross-section 5. Simulated velocities along the right bank (inside of curve, location e) did not match measurements until Cross-section 7. Through the downstream curve, simulated velocities along the right bank (outside of curve, location e) matched measured velocities quite well. Simulated velocities along the left bank (inside of curve, location c) were less than measured velocities, but followed the same decreasing trend with distance through the channel curve. Simulated velocities along the centerline were a bit greater, but within 8 percent of measured velocities. Both the measured and simulated stream-wise velocities along the right bank (outside of curve) increased through the downstream channel curve, while the left bank velocities decreased because of the increasing strength of the secondary circulation through the channel curve.
The simulated stream-wise velocity profiles (using four different mesh sizes) were plotted and compared with the measured velocity profiles from the physical model at Cross-sections 6, 10, and 16 ( Figures A14-A16).
because of the increasing strength of the secondary circulation through the channel curve.
The simulated stream-wise velocity profiles (using four different mesh sizes) were plotted and compared with the measured velocity profiles from the physical model at Cross-sections 6, 10, and 16 ( Figures A14-A16).  Figure A14. Comparisons of simulated and measured stream-wise velocity profiles using fo different mesh sizes at Cross-section 6, located in the upstream channel curve (profile measureme locations "a"-"g" are positioned left to right looking downstream). (g) Figure A14. Comparisons of simulated and measured stream-wise velocity profiles using four different mesh sizes at Cross-section 6, located in the upstream channel curve (profile measurement locations "a"-"g" are positioned left to right looking downstream).  Figure A15. Comparisons of simulated and measured stream-wise velocity profiles using fo different mesh sizes at Cross-section 10, located at the downstream end of the straight transition rea (profile measurement locations "a"-"g" are positioned left to right looking downstream). (g) Figure A15. Comparisons of simulated and measured stream-wise velocity profiles using four different mesh sizes at Cross-section 10, located at the downstream end of the straight transition reach (profile measurement locations "a"-"g" are positioned left to right looking downstream).  Figure A16. Comparisons of simulated and measured stream-wise velocity profiles using fo different mesh sizes at Cross-section 16, located in the downstream channel curve (prof measurement locations "a"-"g" are positioned left to right looking downstream).
In addition to comparisons of the stream-wise profiles, the simulated cross-stream v profiles were also compared with laboratory measurements. The simulated cross-stream vel using the four different mesh sizes, were plotted, and compared with the measured veloc Cross-sections 6, 10, and 16 ( Figures A17-A19). In addition to comparisons of the stream-wise profiles, the simulated cross-stream velocity profiles were also compared with laboratory measurements. The simulated cross-stream velocities, using the four different mesh sizes, were plotted, and compared with the measured velocities at Cross-sections 6, 10, and 16 ( Figures A17-A19). measurement locations "a"-"g" are positioned left to right looking downstream).
In addition to comparisons of the stream-wise profiles, the simulated cross-stream velocity profiles were also compared with laboratory measurements. The simulated cross-stream velocities, using the four different mesh sizes, were plotted, and compared with the measured velocities at Cross-sections 6, 10, and 16 ( Figures A17-A19). The simulated velocity results using the two mesh sizes (492 × 32 × 24 and 492 × 42 × 32) are quite similar. A summary comparison of measured and simulated stream-wise and cross-stream velocities at Cross-section 6 is presented in Table A5 for the finest model mesh size. Similar summary comparisons are presented for Cross-section 10 in Table A6 and for Cross-section 16 in Table A7.  At Cross-section 6, the simulated stream-wise velocity profiles (using the finest 3D mesh) agreed reasonably well with the measurements at locations b, c, d, e, and f ( Figure A14 and Table A5). However, the agreement was not as good at measurement locations a and g on the shallow portions of the left-and right-bank side slopes. The simulated and measured cross-stream velocity profiles were generally in the same direction, but the average of the simulated velocity magnitudes were lower than the measurements ( Figure A17).
At Cross-section 10 (downstream end of the straight transition reach), the simulated and measured stream-wise velocity profiles agreed well at locations b-g ( Figure A15 and Table A6). The patterns of simulated and measured cross-stream velocity profiles match best on the right side of the channel while there is poor agreement on the left side and channel center ( Figure A18). The numerical model predicts that the cross-stream velocities at the left and right banks are toward the channel center, which is consistent with a constricting channel. However, the measured cross-stream velocities are generally toward the left side of the channel, which may be caused by the non-uniform velocity distribution across the upstream boundary of the physical model.
At Cross-section 16 (near the downstream end of the second channel curve), the simulated and measured stream-wise velocity profiles agree well at all measurement locations ( Figure A16 and Table A7). Simulated cross-stream velocity generally did not match measurements ( Figure A19).

Appendix C. Maximum Shear Stress Results
This appendix presents the maximum shear stress results from the 3D numerical modeling and predicted values from the empirical equations (Table A8).