Simulating Laboratory Braided Rivers with Bed-Load Sediment Transport

Numerical models provide considerable assistance in the investigation of complicated processes in natural rivers. In the present study, a physics-based two-dimensional model has been developed to simulate the braiding processes and morphodynamic changes in braided rivers. The model applies the basic hydrodynamic and sediment transport principles with bed morphology deformation and a TVD (Total Variation Diminishing) scheme to predict trans-critical flows and bed morphology deformation. The non-equilibrium transport process of graded bed load sediment is simulated, with non-uniform sediments, secondary flows, and sheltering effects being included. A multiple bed layer technique is adopted to represent the vertical sediment sorting process. The model has been applied to simulate the bed evolution process in an experimental river with bed load transport. Comparisons between the experimental river and predicted river are analysed, including their pattern evolution processes, important braiding phenomena, and statistical characteristics. Avulsion activities have been found in the braiding evolution process, representing the primary ways in which channels form and disappear in braided rivers. The increases in the active braiding intensity and total braiding intensity show similar trends to those observed in the experimental river. Statistical methods are applied to assess the scale-invariant topographic properties of the simulated river and real rivers. The model demonstrated its potential to predict the morphodynamics in natural rivers.


Introduction
Braided rivers are a distinctive river morphology with high flow energy and frequent channel changes.Braiding is one kind of fluvial pattern in these rivers characterized by complex dynamics, in which water and sediment flow are divided into multiple branches, joining and splitting at the nodes and shifting continuously in the floodplain [1].The understanding of morphodynamics in braided rivers is especially important for managing engineering problems such as flood management, riverbank erosion, sedimentation in reservoirs and navigable waterways, and the design and construction of artificial channels and bridges.
Braiding usually occurs with a highly variable discharge, an abundant sediment load, erodible banks, steep valley slopes, and a high width/depth ratio [2].A change in one of these parameters might cause an exchange between meandering and braided patterns such as a large flood, an increased sediment supply, a steeper riverbed, or reduced bank erodibility.Ham and Church studied the lower reach of Chilliwack River in Canada, which is partially braided, and found that a large flood with a return period of about five years can be considered a significant channel-forming event [3].The lower William River in Canada adjusted its channel from a relatively narrow and deep single-channel stream to a thoroughly braided pattern after picking up a 40-fold increase of bed load sediment [4].The Squamish River experienced a downstream sequence of braided-wandering-meandering in a 20 km reach with a slope decrease from 0.0058 to 0.0015 [5].Brierley and Hickin studied flow patterns from minimum flow energy and found that meandering rivers may occur even as flow energy has reached the level for braiding, since scouring has been confined by relatively cohesive or vegetated banks [6].Generally, although one of the factors mentioned above can initiate braided rivers, in nature frequently two or more of them occur simultaneously.Moreover, it is very important to study rivers over an extended time span because their appearance is controlled by the history of changing flow stages, and some geometries arise specially as a result of changing discharge [7].
Numerical models provide useful tools for simulating fluvial processes and morphology in rivers, and, in recent decades, efforts have been made to investigate advanced numerical models to better represent the complex morphodynamic processes in rivers.A few models have been developed to simulate the processes of sediment transport, bed deformation, and sediment downstream fining in rivers (e.g., [8][9][10][11]).Compared with reduced-complexity approaches (e.g., [12,13]), physics-based models provide more detailed process information to better understand natural rivers due to their better representation of hydraulic and morphodynamic processes.Physics-based models have produced braided patterns, bar formation, and avulsion activities observed in natural braided rivers and have made significant progress (e.g., [14,15]).The evolution of wide shallow channels and the behaviour of alternate bars in laboratories have been reproduced by a 2-D model [16,17].Morphodynamic models such as HSTAR (Hydrodynamics and Sediment Transport in Alluvial Rivers) [18], depth-integrated Delft3D [14,19], and 2-D morphodynamic models [20,21] have been applied to simulate large sand-bed braided rivers and have produced morphological units with evolution mechanisms, the characteristic morphology of compound bars, and channels and statistical characteristics that are similar to those of natural sand-bed rivers.However, their development and application are still in an early stage [22], and quantitative analyses and comparisons between simulated and real rivers are still insufficient [14] especially for braided rivers dominated by bed-load sediment transport.
In braided river simulation, previous models have incorporated basic sediment transport theories with essential effects and have shown considerable potential (e.g., [23,24]).The effect of secondary flow has been integrated into the shallow water equations [18,25] or in the sediment transport rate equation [16,17].A TVD (Total Variation Diminishing) scheme, which is a modified Ultimate Quickest scheme, has been employed to solve the flow conservation equations (e.g., [23]), since the scheme is more efficient than the implicit ADI (Alternating Direction Implicit) scheme in predicting trans-critical flows that usually exist in shallow areas of braided rivers.Bank erosion, which is important for the emergence of bars and the initiation of meandering or braiding, has been examined theoretically through the mechanical processes of channels [26,27] and has been proposed to reproduce lateral changes by the existing braided river models (e.g., [17,18]).
However, most of the existing physics-based models for braided river simulation adopted a uniform sediment transport approach (e.g., [14,19,28]), except for the model of Nicholas [18] in which two grain fractions of sand and silt were considered.Nevertheless, non-uniform sediments exist in natural rivers, and different size fractions of sediment particles are transported at different rates [29].Sediment gradation plays an essential role in the bed deformation process, which is especially important for the evolution of local units in natural [30] and laboratory braided rivers [31].In addition, laboratory experiments [32] also show the influence of vertical grain sorting in the active layer, indicating that it is necessary to divide the riverbed into multiple bed layers in river simulation.This mixing layer concept has been adopted by many researchers such as Wu [29] and Van Niekerk et al. [33].Therefore, in order to represent a more realistic morphodynamic process in braided rivers, it is necessary for a numerical model to include the effect of non-uniform sediment transport.
In the current study, a 2-D physics-based model has been developed to simulate an experimental river to study the morphodynamics of braided rivers dominated by bed-load sediment.The model incorporates the basic hydrodynamic and sediment transport theories, considers non-uniform sediments in a fractional way, and adopts multiple bed layers to represent the erosion-deposition process and the sorting effects among fractional sediments in beds.The study aims to investigate the river morphodynamic processes and activities, analyse the statistical characteristics, compare them with the experimental river and natural rivers, and discuss the potential of the model for representing the morphodynamic processes in real rivers.

Numerical Model and Test
Flow, sediment transport, and bed deformation are normally the main processes considered in river morphodynamic simulations.In the present model, flow and sediment transport are calculated by solving the depth-integrated 2D shallow water equations and the non-equilibrium sediment transport equations.

Governing Equations for Flow and Sediment Transport
The flow in rivers is governed by the 2-D shallow water equations, given as follows: ∂p ∂t ∂q ∂t where ξ is the water surface elevation above the datum; p and q are the discharges per unit of width in the x and y directions, respectively; β is the momentum correction factor for a non-uniform vertical velocity profile; U and V are the depth-averaged velocity components in the x and y directions, respectively; f is the Coriolis parameter due to the Earth's rotation; H is the water depth, which is equal to h + ξ, with h as the water depth below the datum; C is the Chézy coefficient, which equals to 18log 10 (12H/k s ), with k s denoting the hydraulic roughness, given according to van Rijn [34]; and ε is the depth-averaged turbulent eddy viscosity.The hydrodynamic model, which is solved by the TVD-MacCormack Scheme, has been validated in previous studies [35,36], showing its efficiency in simulating trans-critical flows.
Assuming that the graded bed material can be divided into N fractions, the transport of the kth fraction can be described by the following equation: where q bk is the bed load transport rate of the kth size fraction; u b is the bed load transport velocity, determined according to the equation of van Rijn [37]; k is the sediment fraction; q b * k is the equilibrium bed load transport rate of the kth size fraction; α bx , α by are the directions of the bed load movement; and L s is the non-equilibrium adaption length, denoting a characteristic distance for the sediment in transport to adjust from a non-equilibrium to an equilibrium state.
The equilibrium bed load transport rate q bk (in m 2 /s) for particles in the range of 0.2 to 2 mm can be calculated by the following equations [38]: where p bk is the percentage of the kth size fraction in the mixing layer of the bed; s is the sediment specific density; and d k is the representative particle diameter of the kth size fraction.The dimensionless transport stage parameter T k and the dimensionless particle parameter D * k are given according to van Rijn [37].

Influence of Bed Slope and Secondary Flow
In natural rivers, when the channel slopes are gentle, the effect of gravity on sediment transport is usually ignored.However, when the riverbed slope is steep, gravity may change the sediment transport capacity and influence bed deformation significantly.According to van Rijn [38], the shear stress with the consideration of the slope effects can be expressed by: where k δ,s is the coefficient of longitudinal slope, with its value equal to sin(θ − δ s )/ sin θ for downslope flow (k δ,s < 1) and sin(θ + δ s )/ sin θ for upslope flow (k δ,s > 1), with θ as the sediment's angle of repose and δ s as the longitudinal slope angle; and k δ,n is the coefficient of the transverse slope, which equals to cos δ n 1 − tan 2 δ n /tan 2 θ 0.5 , with δ n as the transverse slope angle.
The sediment transport rate for a sloping bed is derived as [38]: where α s is the Bagnold slope factor, with its value equal to tan θ/[cos δ s (tan θ ± tan δ s )], where '+' represents upslope flow and '−' represents downslope flow; and q b is the bed load transport rate for a horizontal bed.The bed load transport direction is influenced by the secondary flow and the transverse slope.To consider this effect, the deviation between the bed flow and the sediment transport rate is integrated into the bed load transport rate.Referring to the equations of Jang and Shimizu [16] and van Rijn [38], the sediment transport rate for a transverse slope can be presented as: where, on the right side, the two terms denote the effects of secondary flow and transverse slope, respectively; r s is the radius of curvature of the stream line; N * is the coefficient of the strength of secondary flow, with an value of 7.0 being adopted [39]; and is the mixing coefficient with a value of 1.5 [38].

Areal Fraction and Sheltering Effect for Non-uniform Sediments
To consider the non-uniform sediments in a bed, the forces exerted on graded particles are supposed to be proportional to the projected areas of the particles exposed to flow, and the volumetric percentage in the mixing layer of bed, p bk , can be converted into the areal percentage p ak by the following equation [40]: The sheltering effect in a particle mixture, τ ck , is considered by a correction of the critical shear stress for particle movement initiation of the kth size fraction, giving [41] where d 50 is the median grain size of the sediment mixture.

Bed Deformation and Multiple Bed Layers
During the river evolution process, bed erosion or deposition is determined by the amount of sediment in transport and the sediment transport capacity of the flow.Conversely, the changes of bed morphology and composition will influence the flow velocity and sediment transport rate.The bed deformation can be described as follows: where p is the porosity of the bed material and ∂z b /∂t is the bed elevation change rate.The model supposed that when the bank slope becomes steeper than the sediment's angle of repose, the sediments will immediately slide into neighbouring areas.The amount of sediments deposited on the neighbouring areas is supposed to be equal to what is eroded from the bank.The sediments deposited on the riverbed will be mixed with the upper layer of the bed, and, consequently, the formation and composition of the river bed are renewed.
To consider the influence of bed sediment composition, it is important for numerical models to be capable of resolving the spatial and temporal variation of sediment gradations of the loose layers in the riverbed.This is usually approached by dividing the loose bed into several vertical layers.The present model adopts the three bed layer method of Wei [42], with sediment deposition and erosion calculated in the dynamic layers.More details on the methods of multiple bed layers can be found in Zhou and Lin [43].The sediment transport equations are solved with the Ultimate QUIKEST Scheme, which was developed to simulate 2-D solute transport in coastal and estuarine waters [44].The sediment sorting effect has been investigated [45].More detailed information of the numerical model and model application can be found in Yang [20].

Model Morphodynamic Test
The model for bed load transport was first tested with a laboratory experiment, undertaken by Seal et al. [46] and Paola et al. [47].The experiment was performed in a flume.The test reach was 45 m long and 0.305 m wide, and the initial channel slope was 0.002.A pond located at the downstream end was blocked by a weir, forming a free overfall.The water elevation at the tailgate was set to 0.4 m.The grain size ranged from 0.125 mm to over 64 mm, with two main modes at 0.35 mm and 16 mm.The flow discharge was set to 0.049 m 3 /s, with a sediment feed rate of 11.3 kg/min.The sediment was fed into the channel manually at 1 m downstream of the headgate for 16 h and 50 min.A model was set up with the same boundary condition as the experimental river.The time steps were set to 0.005 s and 0.01 s for the hydrodynamic and sediment computations, respectively.The minimum water depth used to resolve the drying and flooding of the floodplain was 1 mm.
The predicted result together with the experimental data is shown in Figure 1.Generally, the model predicts the bed elevation and the water surface well.As the water and sediment mixture flowed into the channel, sediment particles were deposited onto the bed and caused it to aggrade quickly.The riverbed aggraded vertically with an average speed of 0.0518 m/h, and the downstream head of the newly formed bed moved forward with an average speed of about 2.5 m/h.The aggrading bed was relatively flat with a slope of about 0.2%, which agreed closely with the experimental data.For the longitudinal water surface curve, some oscillations were observed near the upstream inlet when the flow was very shallow, but the oscillations did not develop.The water elevation at 16.8 h was generally simulated well.A hydraulic jump could be observed at the downstream end of the main bed deposit, which was also observed in the experiments.Near the bar, there was a rapid water level change, where the Froude number increased rapidly from 0.2 to 0.9.It illustrates that the hydrodynamic model is capable of coping with sudden hydraulic changes.Figure 2 shows the model-predicted time series of flow velocity, water surface elevation, and bed load transport rate at sequential downstream locations of a 10 m step.It illustrates the variations in local hydraulic conditions as the bed aggrades.Due to a steeper bed slope on the newly formed riverbed, the flow sped up, with a high sediment concentration and a shallower water depth of 0.12 m (compared to the original depth of 0.4 m).Similar changes occurred sequentially along the model river as the newly formed bed grew downstream.

Model Setup
The model was applied to a laboratory experiment conducted by Egozi and Ashmore [31], in which a braided river pattern was generally developed.The river was originally designed to be a generic Froude-scaled model with the prototype being the Sunwapta River in Canada.Froude-scaled models have often been applied in geomorphological studies of gravel braided rivers (e.g., [48][49][50][51]).They satisfy the same Froude number (Fr) as natural rivers yet relax the Reynolds similarity criterion for flows that are hydraulically rough and turbulent.Such a model does not produce exact plan-form similarity for a particular river but yields sampling results applicable to full scale rivers in general [52].
The experiments were performed in an 18 m long and 3 m wide flume, with a constant slope of 0.015.The numerical model was set up by adopting the experimental conditions (Figure 3).The sediments were divided into nine groups, with the size fraction boundaries listed in Table 1.Three bed layers were adopted.Initially the bed was assumed to be flattened, with a trapezoidal straight channel in the middle (Figure 3b).The spatial step was set to 2 cm, while an adaptive time step was used to enhance computational efficiency [35], with the minimum and maximum time steps being 0.001 s and 0.01 s, respectively.The minimum water depth used to resolve the drying and flooding of the floodplain was 1 mm.The model simulation time was approximately 133 h, during which the discharge increased from 1.4 L/s to 2.1 L/s at time t = 71.5 h.The sediments were collected from the outlet and then re-fed to the upstream boundary of inlet, just the same as in the experiment.

Model Setup
The model was applied to a laboratory experiment conducted by Egozi and Ashmore [31], in which a braided river pattern was generally developed.The river was originally designed to be a generic Froude-scaled model with the prototype being the Sunwapta River in Canada.Froude-scaled models have often been applied in geomorphological studies of gravel braided rivers (e.g., [48][49][50][51]).They satisfy the same Froude number (Fr) as natural rivers yet relax the Reynolds similarity criterion for flows that are hydraulically rough and turbulent.Such a model does not produce exact plan-form similarity for a particular river but yields sampling results applicable to full scale rivers in general [52].
The experiments were performed in an 18 m long and 3 m wide flume, with a constant slope of 0.015.The numerical model was set up by adopting the experimental conditions (Figure 3).The sediments were divided into nine groups, with the size fraction boundaries listed in Table 1.Three bed layers were adopted.Initially the bed was assumed to be flattened, with a trapezoidal straight channel in the middle (Figure 3b).The spatial step was set to 2 cm, while an adaptive time step was used to enhance computational efficiency [35], with the minimum and maximum time steps being 0.001 s and 0.01 s, respectively.The minimum water depth used to resolve the drying and flooding of the floodplain was 1 mm.The model simulation time was approximately 133 h, during which the discharge increased from 1.4 L/s to 2.1 L/s at time t = 71.5 h.The sediments were collected from the outlet and then re-fed to the upstream boundary of inlet, just the same as in the experiment.

Model Setup
The model was applied to a laboratory experiment conducted by Egozi and Ashmore [31], in which a braided river pattern was generally developed.The river was originally designed to be a generic Froude-scaled model with the prototype being the Sunwapta River in Canada.Froude-scaled models have often been applied in geomorphological studies of gravel braided rivers (e.g., [48][49][50][51]).They satisfy the same Froude number (Fr) as natural rivers yet relax the Reynolds similarity criterion for flows that are hydraulically rough and turbulent.Such a model does not produce exact plan-form similarity for a particular river but yields sampling results applicable to full scale rivers in general [52].
The experiments were performed in an 18 m long and 3 m wide flume, with a constant slope of 0.015.The numerical model was set up by adopting the experimental conditions (Figure 3).The sediments were divided into nine groups, with the size fraction boundaries listed in Table 1.Three bed layers were adopted.Initially the bed was assumed to be flattened, with a trapezoidal straight channel in the middle (Figure 3b).The spatial step was set to 2 cm, while an adaptive time step was used to enhance computational efficiency [35], with the minimum and maximum time steps being 0.001 s and 0.01 s, respectively.The minimum water depth used to resolve the drying and flooding of the floodplain was 1 mm.The model simulation time was approximately 133 h, during which the discharge increased from 1.4 L/s to 2.1 L/s at time t = 71.5 h.The sediments were collected from the outlet and then re-fed to the upstream boundary of inlet, just the same as in the experiment.

River Evolution Processes and Phenomena
The time serial pictures denoting the evolution process of the model-predicted river are shown in Figure 4.It can be observed that the development of a braided pattern in the predicted river shares many common properties with the experimental river (Figure 5, activity number identified with Figure 4).Most of the morphologic elements and activities observed in the experimental river [31] have been found in the predicted river.
In the first few hours of model simulation, alternate bars and channels were formed in the upstream reach, mainly distributed in the initially cut channel (hour 8).Next, one main channel began to take shape under the strengthening connection of the alternate channels by incision, bank erosion, and channel migration (hour 12).In this channel, the water depth increased further along with increasing sinuosity (hour 19, point 1).Intense erosion occurred along the bank of the initially cut channel, and the channel migrated outward.As the channel became more sinuous, the water overflowed from the banks of certain local bends and spread across the unchannelled margins downstream of the bends (hour 24, point 2), which were called 'pseudoanabranches' and sometimes grew into new channels.At certain local bends, the channel was dissected by the upstream tributary, and one bifurcation was formed (hour 24, point 3).At the bifurcation of the newly avulsed areas, the flow rate decreased in the original main channel and increased in the new channel, and the new channel consequently became more dominant, whilst the old channel died out gradually.This activity can be called avulsion by incision [15], through which the main channel tends to find a straighter pathway by incising a new course.
The braiding configuration was developed by hour 27.Some active channels reconnected with each other again by overflow or new channel erosion (hour 27, point 5), while some secondary channels were gradually abandoned (hour 27, point 4).As some channels became more sinuous, water overflowed at the bends, forming new channels (hour 24, point 6).Flow concentrated in the main channel, and accordingly it formed another channel bend (hour 36, point 7).Certain channels were choked by sediment deposition (hour 38, point 8), which can be classified as avulsion by progradation [15].By hour 40, the upstream reach had widened significantly by extending downstream and into the original margin areas, evolving into a new stage.One of the bifurcated  The time serial pictures denoting the evolution process of the model-predicted river are shown in Figure 4.It can be observed that the development of a braided pattern in the predicted river shares many common properties with the experimental river (Figure 5, activity number identified with Figure 4).Most of the morphologic elements and activities observed in the experimental river [31] have been found in the predicted river.
In the first few hours of model simulation, alternate bars and channels were formed in the upstream reach, mainly distributed in the initially cut channel (hour 8).Next, one main channel began to take shape under the strengthening connection of the alternate channels by incision, bank erosion, and channel migration (hour 12).In this channel, the water depth increased further along with increasing sinuosity (hour 19, point 1).Intense erosion occurred along the bank of the initially cut channel, and the channel migrated outward.As the channel became more sinuous, the water overflowed from the banks of certain local bends and spread across the unchannelled margins downstream of the bends (hour 24, point 2), which were called 'pseudoanabranches' and sometimes grew into new channels.At certain local bends, the channel was dissected by the upstream tributary, and one bifurcation was formed (hour 24, point 3).At the bifurcation of the newly avulsed areas, the flow rate decreased in the original main channel and increased in the new channel, and the new channel consequently became more dominant, whilst the old channel died out gradually.This activity can be called avulsion by incision [15], through which the main channel tends to find a straighter pathway by incising a new course.
The braiding configuration was developed by hour 27.Some active channels reconnected with each other again by overflow or new channel erosion (hour 27, point 5), while some secondary channels were gradually abandoned (hour 27, point 4).As some channels became more sinuous, water overflowed at the bends, forming new channels (hour 24, point 6).Flow concentrated in the main channel, and accordingly it formed another channel bend (hour 36, point 7).Certain channels were choked by sediment deposition (hour 38, point 8), which can be classified as avulsion by progradation [15].By hour 40, the upstream reach had widened significantly by extending downstream and into the original margin areas, evolving into a new stage.One of the bifurcated channels captured a large amount of flow from another and finally evolved into a deep channel, followed by the disappearance of the original dominant channel of a bifurcation (hour 44, point 9).The main channel shifted to one of its tributaries as the upstream flow conditions changed (hour 44, point 10).Certain previous abandoned channels (hour 38, point 8) were regenerated by channel annexation (hour 47, point 11).In the predicted river, a cycle of straight channels transitioning to sinuous channels and then to avulsions is an important morphodynamic process during the evolutionary process of the braided river.These activities play an important role in maintaining the dynamic braided form.

Properties of a Typical Braided River
The model riverbed was initially flat with a straight channel in the middle of the flume, and it gradually evolved into a braided river.A typical braided river was generated after 41 h, with the distributions of water depth, Froude number, bed shear stress, bed load concentration, and sediment distribution shown in Figure 6.Key channel nodes, including bars, confluences, bifurcations, and a series of sequential pool-bar units with repeated division and joining of channels can be observed (Figure 6a).These features share many similarities with the experimental river and natural rivers (e.g., [7,31,53]).
For most of the active channels, the Froude number was between 0.6 and 0.7 (Figure 6b).In some local areas, especially in the areas where confluence and bifurcation occurred, the Froude number varied in a range of 0.7 to 1.0.The distribution of bed shear stress has been shown to be relative to bed grain size distribution and water depth (Figure 6c).Higher shear stress values occurred mainly in the main channel and bifurcation areas, and at times they also existed along channel edges.
An active channel is defined as a channel with visible bed material movement [54].In the simulated river, if the bed load concentration is higher than a certain value (e.g., 6 g/m 3 ), it can be considered 'active'; otherwise, it is considered 'non-active'.In that sense, a main active channel

Properties of a Typical Braided River
The model riverbed was initially flat with a straight channel in the middle of the flume, and it gradually evolved into a braided river.A typical braided river was generated after 41 h, with the distributions of water depth, Froude number, bed shear stress, bed load concentration, and sediment distribution shown in Figure 6.Key channel nodes, including bars, confluences, bifurcations, and a series of sequential pool-bar units with repeated division and joining of channels can be observed (Figure 6a).These features share many similarities with the experimental river and natural rivers (e.g., [7,31,53]).
For most of the active channels, the Froude number was between 0.6 and 0.7 (Figure 6b).In some local areas, especially in the areas where confluence and bifurcation occurred, the Froude number varied in a range of 0.7 to 1.0.The distribution of bed shear stress has been shown to be relative to bed grain size distribution and water depth (Figure 6c).Higher shear stress values occurred mainly in the main channel and bifurcation areas, and at times they also existed along channel edges.
An active channel is defined as a channel with visible bed material movement [54].In the simulated river, if the bed load concentration is higher than a certain value (e.g., 6 g/m 3 ), it can be considered 'active'; otherwise, it is considered 'non-active'.In that sense, a main active channel formed in the river with the highest sediment concentration, and other active channels joined and bifurcated from it (Figure 6d), which was similar to the experimental river.Due to the special sediment load method in the experiment in which 'the sand leaving the end of the flume was recirculated from the tail tank and fed back into the flume' [31], finer sediments were more involved in the circulation because it was easier for them to be transported; therefore, the sediments were often coarser on the two sides of a channel than in the thalweg (Figure 6e).
formed in the river with the highest sediment concentration, and other active channels joined and bifurcated from it (Figure 6d), which was similar to the experimental river.Due to the special sediment load method in the experiment in which 'the sand leaving the end of the flume was recirculated from the tail tank and fed back into the flume' [31], finer sediments were more involved in the circulation because it was easier for them to be transported; therefore, the sediments were often coarser on the two sides of a channel than in the thalweg (Figure 6e).

Sensitivity Analysis
A sensitivity analysis was performed to investigate the effects of secondary flow and grid size on the model performance.Without considering the effect of secondary flow (Figure 7a), a river with alternating pool-bar units also formed, and a braided pattern was generated at hour 41.This river shared a similar channel pattern with the river considering the effect of secondary flow (Figure 7b), yet it generated fewer but wider channels with a slower development rate, illustrating the role that the secondary flow plays in transverse sediment transport and channel migration.
Figure 8 shows the braiding configurations with grids of (a) 0.03 × 0.03 m 2 , (b) 0.025 × 0.025 m 2 , and (c) 0.02 × 0.02 m 2 at hour 27.Compared with the river with a cell size of 0.02 × 0.02 m 2 , the river with a lower resolution of 0.03 × 0.03 m 2 produced fewer and wider channels.In contrast, the river with a middle cell resolution of 0.025 × 0.025 m 2 generated a similar braided pattern to the river with a cell size of 0.02 × 0.02 m 2 .The channel numbers and widths of these two rivers were comparable with each other.

Sensitivity Analysis
A sensitivity analysis was performed to investigate the effects of secondary flow and grid size on the model performance.Without considering the effect of secondary flow (Figure 7a), a river with alternating pool-bar units also formed, and a braided pattern was generated at hour 41.This river shared a similar channel pattern with the river considering the effect of secondary flow (Figure 7b), yet it generated fewer but wider channels with a slower development rate, illustrating the role that the secondary flow plays in transverse sediment transport and channel migration.
Figure 8 shows the braiding configurations with grids of (a) 0.03 × 0.03 m 2 , (b) 0.025 × 0.025 m 2 , and (c) 0.02 × 0.02 m 2 at hour 27.Compared with the river with a cell size of 0.02 × 0.02 m 2 , the river with a lower resolution of 0.03 × 0.03 m 2 produced fewer and wider channels.In contrast, the river with a middle cell resolution of 0.025 × 0.025 m 2 generated a similar braided pattern to the river with a cell size of 0.02 × 0.02 m 2 .The channel numbers and widths of these two rivers were comparable with each other.

Statistical Characteristics
Due to the nature of braided rivers with frequent channel migrations, statistical methods are usually employed to assess the characteristics of braided rivers [55].Given that braided rivers in nature exhibit some intrinsic characteristics common to all braided rivers, several methods have been proposed to evaluate the simulation capability of physical and numerical models (e.g., [55,56]).In the present study, three statistical methods, including braiding intensity, state-space plots, and transect topography and slope frequency, were applied to evaluate the model-predicted river.

Braiding Intensity
Braiding intensity has been proposed to measure the channel intensity in braided rivers [54,57,58].Total braiding intensity (BIT) (defined as the number of total wetted channels counted and averaged over a number of cross sections of a river), active braiding intensity (BIA) (defined as the number of channels with visible bed load sediment transport counted and averaged over a number of cross sections of a river), and their ratio (BIT/BIA) were calculated (Figure 9) and compared with the experimental results (Table 2).In general, the model successfully predicted the response of the

Statistical Characteristics
Due to the nature of braided rivers with frequent channel migrations, statistical methods are usually employed to assess the characteristics of braided rivers [55].Given that braided rivers in nature exhibit some intrinsic characteristics common to all braided rivers, several methods have been proposed to evaluate the simulation capability of physical and numerical models (e.g., [55,56]).In the present study, three statistical methods, including braiding intensity, state-space plots, and transect topography and slope frequency, were applied to evaluate the model-predicted river.

Braiding Intensity
Braiding intensity has been proposed to measure the channel intensity in braided rivers [54,57,58].Total braiding intensity (BIT) (defined as the number of total wetted channels counted and averaged over a number of cross sections of a river), active braiding intensity (BIA) (defined as the number of channels with visible bed load sediment transport counted and averaged over a number of cross sections of a river), and their ratio (BIT/BIA) were calculated (Figure 9) and compared with the experimental results (Table 2).In general, the model successfully predicted the response of the

Statistical Characteristics
Due to the nature of braided rivers with frequent channel migrations, statistical methods are usually employed to assess the characteristics of braided rivers [55].Given that braided rivers in nature exhibit some intrinsic characteristics common to all braided rivers, several methods have been proposed to evaluate the simulation capability of physical and numerical models (e.g., [55,56]).In the present study, three statistical methods, including braiding intensity, state-space plots, and transect topography and slope frequency, were applied to evaluate the model-predicted river.

Braiding Intensity
Braiding intensity has been proposed to measure the channel intensity in braided rivers [54,57,58].Total braiding intensity (BI T ) (defined as the number of total wetted channels counted and averaged over a number of cross sections of a river), active braiding intensity (BI A ) (defined as the number of channels with visible bed load sediment transport counted and averaged over a number of cross sections of a river), and their ratio (BI T /BI A ) were calculated (Figure 9) and compared with the experimental results (Table 2).In general, the model successfully predicted the response of the channel pattern to increasing discharge, with the total and active braiding intensities increased to new stages.The maximum BI T values were greater than twice the maximum BI A values, which were similar to those observed in the experiments [31].Compared with the experimental river, the model-predicted river developed nearly the same BI A and a slightly lower BI T in Experiment 7.However, in Experiment 8, the differences between the model-predicted and the measured BI A and BI T values increased after hour 100.This is because several channels co-existed and the difference between the main channel and other channels became less obvious.
Water 2017, 9, 686 12 of 17 channel pattern to increasing discharge, with the total and active braiding intensities increased to new stages.The maximum BIT values were greater than twice the maximum BIA values, which were similar to those observed in the experiments [31].Compared with the experimental river, the modelpredicted river developed nearly the same BIA and a slightly lower BIT in Experiment 7.However, in Experiment 8, the differences between the model-predicted and the measured BIA and BIT values increased after hour 100.This is because several channels co-existed and the difference between the main channel and other channels became less obvious.In both flow stages, the active braiding intensity, BIA, developed quickly to a stable value (Figure 9).In the first stage, it took slightly longer to stabilize, but, in the second stage, it took only a couple of hours.In contrast, BIT took a much longer time to increase progressively to a stable state (approximately 40 h).The quick stabilization of BIA was explained by Egozi and Ashmore [31] as follows: the increase in BIA requires only a sufficiently large flow within an existing channel to mobilize the bed material, and this can be accomplished very quickly; however, the increase in BIT needs time for the river to erode new channels, resulting in a much longer time for its stabilization.The longer duration for the stabilization of BIA and BIT during the first flow stage resulted from the fact that it took some time to develop a braided pattern from the initially straight channel.The variation in BIA/BIT shows a consistent trend progressively approaching nearly 0.4 (Figure 9), which shared similar stable values of BIA / BIT with the experimental river but with a slightly narrower range (Table 2).

State Space Plots
Natural braided rivers have shown an intertwining effect of braiding, represented by the repeated division and joining of channels with bars between them [7].Flow at one location affects the downstream morphology, and braided rivers show recurring geometric and topographic  In both flow stages, the active braiding intensity, BI A , developed quickly to a stable value (Figure 9).In the first stage, it took slightly longer to stabilize, but, in the second stage, it took only a couple of hours.In contrast, BI T took a much longer time to increase progressively to a stable state (approximately 40 h).The quick stabilization of BI A was explained by Egozi and Ashmore [31] as follows: the increase in BI A requires only a sufficiently large flow within an existing channel to mobilize the bed material, and this can be accomplished very quickly; however, the increase in BI T needs time for the river to erode new channels, resulting in a much longer time for its stabilization.The longer duration for the stabilization of BI A and BI T during the first flow stage resulted from the fact that it took some time to develop a braided pattern from the initially straight channel.The variation in BI A /BI T shows a consistent trend progressively approaching nearly 0.4 (Figure 9), which shared similar stable values of BI A / BI T with the experimental river but with a slightly narrower range (Table 2).

State Space Plots
Natural braided rivers have shown an intertwining effect of braiding, represented by the repeated division and joining of channels with bars between them [7].Flow at one location affects the downstream morphology, and braided rivers show recurring geometric and topographic characteristics.A method named 'state-space plot' has been employed to evaluate the spatial recurring patterns of braided rivers (e.g., [56]).The state-space plot is constructed using a series of values of one essential variable by plotting each value versus the previous value at its neighbouring upstream cross section, and when it moves from one point to the successive points and eventually forms a closed loop, the information on spatial order is recorded.The total channel widths, sequential scour depths at channel confluences, and heights of successive bars have been proposed to analyse the spatial repeating pattern of braided rivers [21,55,59,60].
In the present study, the maximum scour depth in each cross-section was chosen for state-space plotting.With the data of a fully developed braided river, state-space plots were constructed by plotting the scour depth versus the preceding scour depth, both normalised to the overall scour depth.Figure 10 shows the state space plots of the predicted river and a laboratory river [60].Clusters of loops can be seen in the plots of the predicted river and the laboratory river, indicating that they share a common statistically spatial repeating pattern.However, the predicted river has a thinner pattern than the laboratory river, indicating that, in the laboratory river, the channel scour depth changes in a more abrupt way than in the predicted river.Nonetheless, crossing lines were more common for the laboratory river than for the predicted river.They were considered to originate from some statistic behaviour of the rivers [59].
Water 2017, 9, 686 13 of 17 characteristics.A method named 'state-space plot' has been employed to evaluate the spatial recurring patterns of braided rivers (e.g., [56]).The state-space plot is constructed using a series of values of one essential variable by plotting each value versus the previous value at its neighbouring upstream cross section, and when it moves from one point to the successive points and eventually forms a closed loop, the information on spatial order is recorded.The total channel widths, sequential scour depths at channel confluences, and heights of successive bars have been proposed to analyse the spatial repeating pattern of braided rivers [21,55,59,60].
In the present study, the maximum scour depth in each cross-section was chosen for state-space plotting.With the data of a fully developed braided river, state-space plots were constructed by plotting the scour depth versus the preceding scour depth, both normalised to the overall scour depth.Figure 10 shows the state space plots of the predicted river and a laboratory river [60].Clusters of loops can be seen in the plots of the predicted river and the laboratory river, indicating that they share a common statistically spatial repeating pattern.However, the predicted river has a thinner pattern than the laboratory river, indicating that, in the laboratory river, the channel scour depth changes in a more abrupt way than in the predicted river.Nonetheless, crossing lines were more common for the laboratory river than for the predicted river.They were considered to originate from some statistic behaviour of the rivers [59].

Transect Topography and Slope Frequency
The riverbed topography structure can be described by cross-section surveys such as the transect topography and the cumulative frequency of lateral slopes [55].The former provides direct transect distributions of bed topographic highs and lows.The latter estimates the frequency of bed lateral slopes relative to a reference level; for example, the median elevation.They can be employed to assess the similarities between the simulated and natural rivers.
The predicted river with a fully developed braided pattern was analysed to determine the structure of the transect topography.The results of two representative cross-sections are shown in Figure 11, together with those from the Sunwapta River and a laboratory river [55], with the Sunwapta River being the prototype of the experimental river simulated in the present study.It can be seen that, in the cross-sections shown in Figure 11, the channels are relatively narrow and deep, while the bars are relatively wide and gentle, with the maximum scour depth of the channels exceeding the maximum deposition height of the bars.These are common characteristics for the model-predicted river and the prototype natural and laboratory rivers.However, compared to the Sunwapta River and laboratory river, the predicted river presents a smoother cross-sectional bed elevation.This indicates that the riverbed of the simulated river is more ideal than those of real rivers.

Transect Topography and Slope Frequency
The riverbed topography structure can be described by cross-section surveys such as the transect topography and the cumulative frequency of lateral slopes [55].The former provides direct transect distributions of bed topographic highs and lows.The latter estimates the frequency of bed lateral slopes relative to a reference level; for example, the median elevation.They can be employed to assess the similarities between the simulated and natural rivers.
The predicted river with a fully developed braided pattern was analysed to determine the structure of the transect topography.The results of two representative cross-sections are shown in Figure 11, together with those from the Sunwapta River and a laboratory river [55], with the Sunwapta River being the prototype of the experimental river simulated in the present study.It can be seen that, in the cross-sections shown in Figure 11, the channels are relatively narrow and deep, while the bars are relatively wide and gentle, with the maximum scour depth of the channels exceeding the maximum deposition height of the bars.These are common characteristics for the model-predicted river and the prototype natural and laboratory rivers.However, compared to the Sunwapta River and laboratory river, the predicted river presents a smoother cross-sectional bed elevation.This indicates that the riverbed of the simulated river is more ideal than those of real rivers.The slope frequency distribution curves were constructed based on the cross-section slopes in the model predicted reach of 4 to 17 m, containing the data of 650 cross-sections.Figure 12 shows the results of the predicted river at hour 41, the Sunwapta River, and a laboratory river [55], all of which present similar varying trends.The curves corresponding to the channels are significantly wider than the bars, illustrating that the bars exhibit statistically gentler elevation changes than do the channels.This is consistent with the findings from the bed-transect topographic distributions (Figure 11).The steep sides in the channels were produced by fast flow and energetic erosion activity, while an opposite situation occurred on the bars.In addition, the wider slope frequency domains of both the channels and the bars for the predicted river indicate that statistically the predicted river channels are slightly steeper than those of the Sunwapta River and the laboratory river.

Discussion and Conclusions
A physics-based numerical model was developed to simulate the morphodynamic processes in laboratory braided rivers with non-uniform bed load transport.The model applied the basic hydrodynamic and sediment transport principles with a fractional method and multiple bed layers.It employed the TVD-MacCormack Scheme for solving the hydrodynamic equations of trans-critical flows, and was tested with a sediment aggradation case and found to be capable of effectively The slope frequency distribution curves were constructed based on the cross-section slopes in the model predicted reach of 4 to 17 m, containing the data of 650 cross-sections.Figure 12 shows the results of the predicted river at hour 41, the Sunwapta River, and a laboratory river [55], all of which present similar varying trends.The curves corresponding to the channels are significantly wider than the bars, illustrating that the bars exhibit statistically gentler elevation changes than do the channels.This is consistent with the findings from the bed-transect topographic distributions (Figure 11).The steep sides in the channels were produced by fast flow and energetic erosion activity, while an opposite situation occurred on the bars.In addition, the wider slope frequency domains of both the channels and the bars for the predicted river indicate that statistically the predicted river channels are slightly steeper than those of the Sunwapta River and the laboratory river.The slope frequency distribution curves were constructed based on the cross-section slopes in the model predicted reach of 4 to 17 m, containing the data of 650 cross-sections.Figure 12 shows the results of the predicted river at hour 41, the Sunwapta River, and a laboratory river [55], all of which present similar varying trends.The curves corresponding to the channels are significantly wider than the bars, illustrating that the bars exhibit statistically gentler elevation changes than do the channels.This is consistent with the findings from the bed-transect topographic distributions (Figure 11).The steep sides in the channels were produced by fast flow and energetic erosion activity, while an opposite situation occurred on the bars.In addition, the wider slope frequency domains of both the channels and the bars for the predicted river indicate that statistically the predicted river channels are slightly steeper than those of the Sunwapta River and the laboratory river.

Discussion and Conclusions
A physics-based numerical model was developed to simulate the morphodynamic processes in laboratory braided rivers with non-uniform bed load transport.The model applied the basic hydrodynamic and sediment transport principles with a fractional method and multiple bed layers.It employed the TVD-MacCormack Scheme for solving the hydrodynamic equations of trans-critical flows, and was tested with a sediment aggradation case and found to be capable of effectively

Discussion and Conclusions
A physics-based numerical model was developed to simulate the morphodynamic processes in laboratory braided rivers with non-uniform bed load transport.The model applied the basic hydrodynamic and sediment transport principles with a fractional method and multiple bed layers.It employed the TVD-MacCormack Scheme for solving the hydrodynamic equations of trans-critical flows, and was tested with a sediment aggradation case and found to be capable of effectively predicting the flow and bed deformation.
With the same boundary conditions as an experimental river, the model successfully reproduced the evolution processes and phenomena of braiding, avulsion activities, and the responses of the channel pattern to increased discharge.The incorporation of the effect of secondary flow played an essential role in producing a more realistic braided pattern.Certain geomorphologic phenomena and processes observed in the experimental river and natural rivers have been found in the evolutionary process of the predicted river.A cycle of straight channels transitioning to sinuous channels and then to avulsions was an important morphodynamic process for the river to maintain a dynamic braided form.Avulsion activities with similar evolution processes to natural rivers have been found in both the experimental and predicted rivers.A high shear stress usually occurred in areas where the flow was energetic, erosion was possible, and the riverbed sediment was mainly composed of coarse bed material.The changes in the braiding intensities illustrated trends similar to those in the experimental river in the process of braiding generation and in response to increased discharge.Active braiding intensity developed quickly to a stable value, whereas total braiding intensity increased to a stable state gradually due to the time consumed by cutting new channels.The predicted river shared common statistical topographic characteristics with natural and laboratory rivers, showing a recurring spatial pattern of shallow and deep topography and statistical transect gentle bars and steep channels.
Generally, the results presented above demonstrate the potential for the present model to represent the morphodynamic processes of natural braided rivers.Nevertheless, it is necessary to gradually improve the understanding and capability of the model in conjunction with field research and laboratory experiments.
0.4 m).Similar changes occurred sequentially along the model river as the newly formed bed grew downstream.

Figure 1 .
Figure 1.Predicted bed elevations and water surface of the experiment of Seal et al. [46].

Figure 1 .
Figure 1.Predicted bed elevations and water surface of the experiment of Seal et al. [46].

Figure 1 .
Figure 1.Predicted bed elevations and water surface of the experiment of Seal et al. [46].
large amount of flow from another and finally evolved into a deep channel, followed by the disappearance of the original dominant channel of a bifurcation (hour 44, point 9).The main channel shifted to one of its tributaries as the upstream flow conditions changed (hour 44, point 10).Certain previous abandoned channels (hour 38, point 8) were regenerated by channel annexation (hour 47, point 11).In the predicted river, a cycle of straight channels transitioning to sinuous channels and then to avulsions is an important morphodynamic process during the evolutionary process of the braided river.These activities play an important role in maintaining the dynamic braided form.

Figure 4 .
Figure 4. Morphologic changes and activities during the evolution process of the model predicted river (water depth/mm).1: increasing sinuosity of the main channel; 2: the generation and obliteration of pseudo-anabranches; 3: avulsion by incision; 4: abandonment of the secondary channel; 5: reconnecting of two active channels; 6: avulsion by overflow; 7: concentration of the flow in the main channel; 8: avulsion by progradation; 9: water diversion from one secondary channel to an existing or new channel; 10: shifting of the main channel; and 11: regeneration of an abandoned channel.

Figure 4 .
Figure 4. Morphologic changes and activities during the evolution process of the model predicted river (water depth/mm).1: increasing sinuosity of the main channel; 2: the generation and obliteration of pseudo-anabranches; 3: avulsion by incision; 4: abandonment of the secondary channel; 5: reconnecting of two active channels; 6: avulsion by overflow; 7: concentration of the flow in the main channel; 8: avulsion by progradation; 9: water diversion from one secondary channel to an existing or new channel; 10: shifting of the main channel; and 11: regeneration of an abandoned channel.

Figure 5 .
Figure 5. Morphologic changes and activities during the evolution process of the experimental river with the activity numbers identical to those in Figure 4.

Figure 5 .
Figure 5. Morphologic changes and activities during the evolution process of the experimental river with the activity numbers identical to those in Figure 4.

Figure 7 .
Figure 7. Braiding configurations (a) without and (b) with the effect of secondary flow.

Figure 7 . 17 Figure 7 .
Figure 7. Braiding configurations (a) without and (b) with the effect of secondary flow.

Figure 9 .
Figure 9. Variations in total braiding intensity (BI T ), active braiding intensity (BI A ), and BI A /BI T during sequential flow stages.

Figure 10 .
Figure 10.State-space plots for (a) the model-predicted river and (b) a laboratory river.

Figure 10 .
Figure 10.State-space plots for (a) the model-predicted river and (b) a laboratory river.

Water 2017, 9 , 686 14 of 17 Figure 11 .
Figure 11.Deviations from the elevation median of cross-sections for (a) the model predicted river; (b) the Sunwapta River; and (c) a laboratory river.The dotted line represents the elevation median.

Figure 12 .
Figure 12.Cumulative distributions of the lateral slopes of the bars and channels for (a) the model predicted river; (b) the Sunwapta River; and (c) a laboratory river.

Figure 11 .
Figure 11.Deviations from the elevation median of cross-sections for (a) the model predicted river; (b) the Sunwapta River; and (c) a laboratory river.The dotted line represents the elevation median.

Water 2017, 9 , 686 14 of 17 Figure 11 .
Figure 11.Deviations from the elevation median of cross-sections for (a) the model predicted river; (b) the Sunwapta River; and (c) a laboratory river.The dotted line represents the elevation median.

Figure 12 .
Figure 12.Cumulative distributions of the lateral slopes of the bars and channels for (a) the model predicted river; (b) the Sunwapta River; and (c) a laboratory river.

Figure 12 .
Figure 12.Cumulative distributions of the lateral slopes of the bars and channels for (a) the model predicted river; (b) the Sunwapta River; and (c) a laboratory river.

Table 1 .
Sediment fractions in the experimental river and the numerical model.

Table 1 .
Sediment fractions in the experimental river and the numerical model.

Table 2 .
The braiding intensities of the model predicted river and the experimental river.