## 2 Materials and Methods
Because it is not possible to compute the flow dynamics about each individual grain of sediment, an empirical model must be used. The model used here is based on Mastbergen and Van den Berg. The first step to computing the critical Shields parameter is calculating the dimensionless parameter \(d_{*;i}\):  
\[d_{*;i}=d_{i}\left[\frac{\rho_{f}\left(\rho_{i}-\rho_{f}\right)\|g\|}{\mu_{f}^ {2}}\right]^{1/3} \tag{7}\]  
where \(\rho_{i}\) is the density of the sediment species I, \(\rho_{f}\) is the fluid density, \(d_{i}\) is the diameter, \(\mu_{f}\) is the dynamic viscosity of fluid, and \(\|g\|\) is the magnitude if the acceleration of gravity g.  
From this, the dimensionless critical Shields parameter is computed using the Soulsby-Whitehouse equation:  
\[\theta_{cr,i}=\frac{0.3}{1+1.2d_{*;i}}+0.055[1-\exp(-0.02d_{*;i})] \tag{8}\]  
The critical Shields parameter can be modified for sloping surfaces to include the angle of repose. The modification further alters \(\theta_{cr,i}\)  
\[\theta_{cr,i}^{\prime}=\theta_{cr,i}\frac{\cos\psi\sin\beta+\sqrt{\cos^{2} \beta\tan^{2}\ \phi_{i}-\sin^{2}\psi\sin^{2}\beta}}{\tan\ \phi_{i}} \tag{9}\]  
where \(\beta\) is the angle of slope of bed, \(\varphi\) is the user-defined angle of repose for sediment species i (default is 32), and \(\psi\) is the angle between the flow and the upslope direction [17]. The local Shields parameter is computed based on the local bed shear stress, \(\tau\):  
\[\theta_{i}=\frac{\tau}{\|g\|_{d_{i}}\left(\rho_{i}-\rho_{f}\right)} \tag{10}\]  
where \(\tau\) is calculated using the law of the wall for 3D turbulent flow and the quadratic law of bottom shear stress for shallow-water turbulent flow, respectively, with consideration of the bed surface roughness. The settling velocity equation proposed by [PERSON] [19] is used:  
\[\frac{\partial U_{i}}{\partial t}+\frac{1}{V_{f}}\left(U_{i}A_{j}\frac{ \partial u}{\partial X_{j}}\right)=-\frac{1}{\rho}\frac{\partial P}{\partial X _{i}}+G_{i}+f_{i} \tag{2}\]0  
Bed-load transport is a mode of sediment transport caused by rolling or bouncing over the surface of the packed bed of sediment. The volumetric transport rate of sediment per width of bed is as follows [20]:  
\[\frac{\partial U_{i}}{\partial t}+\frac{1}{V_{f}}\left(U_{i}A_{j}\frac{ \partial u}{\partial X_{j}}\right)=-\frac{1}{\rho}\frac{\partial P}{\partial X _{i}}+G_{i}+f_{i} \tag{2}\]1  
where the \(\beta_{VT;i}\) coefficient is 0.053. \(c_{b,i}\) is the volume fraction of species i in the bed material. For each species, the suspended sediment concentration is calculated by solving its own transport equation,  
\[\frac{\partial U_{i}}{\partial t}+\frac{1}{V_{f}}\left(U_{i}A_{j}\frac{ \partial u}{\partial X_{j}}\right)=-\frac{1}{\rho}\frac{\partial P}{\partial X _{i}}+G_{i}+f_{i} \tag{2}\]2Here, \(C_{s,i}\) is the suspended sediment mass concentration of species i, which is defined as the sediment mass per volume of fluid-sediment mixture; \(D\) is the diffusivity; and \(u_{s,i}\) is the suspended sediment velocity. The FAVOR (Fractional Area/Volume Obstacle Representation) method for the description of the geometry and the Volume-of-Fluid (VOF) method as originally described by [PERSON] [21] for tracking fluid interfaces are used in the model. Typically, the fractional volume is represented by the quantity. In interior regions of liquid, the value of \(F\) would be 1.0, while outside of the liquid, in regions of gas (air, for example), the value of \(F\) is zero. The location of a free surface is where \(F\) changes from 0.0 to 1.0. Thus, any element having an \(F\) value lying between 0.0 and 1.0 must contain a surface. In addition, the normal to the surface can be calculated from the direction in which \(F\) changes most rapidly, applying boundary conditions to the surface.